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SECTION  1 

INTRODUCTION  AND  SUMMARY 


This  report  documents  a computer  code  for  dynamic  analyses  of  impact  and 
explosive  detonation  problems.  The  code,  EPIC-2  (Elastic -Plastic  hnpact 
Computations  in  2_  dimens  ions)  is  applicable  for  axi-symmetric  and  plane 
strain  problems.  It  also  has  the  ability  to  handle  the  effect  of  spin  for 
the  axi-symmetric  case.  It  is  based  on  a Lagrangian  finite  element  formu- 
lation where  the  equations  of  motion  are  integrated  directly,  rather  than 
4 through  the  traditional  stiffness  matrix  approach.  Nonlinear  material 

strength  and  compressibility  effects  are  included  to  account  for  elastic- 
plastic  flow  and  wave  propagation.  An  option  to  include  high  explosives  is 
also  available.  Although  the  code  is  arranged  to  provide  solutions  for  pro- 
jectile-target impact,  and  explosive  detonation  problems,  the  basic  formula- 
tion is  valid  for  a wide  range  of  problems  involving  dynamic  responses  of 
continuous  media, 
i 

The  EPIC-2  code  has  material  descriptions  which  include  strain  hardening, 
strain  rate  effects,  thermal  softening  and  fracture.  Geometry  generators 
are  included  to  quickly  generate  flat  plates,  spheres  and  rods  with  blunt, 
rounded  or  conical  nose  shapes.  It  has  the  capacity  to  include  multiple  slid- 
ing surfaces,  it  is  restartable,  and  it  provides  plots  of  initial  and  deformed 
geometry  as  well  as  strain,  stress,  pressure  and  velocity  fields.  Time- 
dependent  plots  of  various  system  parameters  are  also  available. 

A desirable  characteristic  of  this  technique  is  that  there  is  no  need  to  pro- 
vide an  orderly  arrangement  of  nodes  as  is  required  for  finite-difference 
techniques.  Complex  geometrical  shapes  can  be  represented  simply  by 
providing  an  adequate  assemblage  of  elements  to  represent  the  geometry. 
Various  boundary  conditions  can  also  be  represented  in  a straightforward 
manner.  Another  desirable  characteristic  is  that  this  technique  is  formulated 
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1 


j 

for  a triangular  element,  which  is  well  suited  to  represent  the  severe  dis- 
tortions which  often  occur  during  high-velocity  impact.  A triangular  element 
also  provides  a state  of  constant  strain  such  that  all  material  in  an  element 
behaves  uniformly.  Hus  allows  for  an  accurate  and  convenient  selection  of 
a constant  stress  within  the  element. 

TTiis  report  includes  a complete  description  and  formulation  of  the  EPIC-2 
computer  code.  Detailed  instructions  for  using  this  code  are  also  provided. 
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SECTION  2 
FORMULATION 


The  EPIC-2  computational  technique  for  axi-symmetric  impact  problems  is 
shown  in  Figure  1.  The  first  step  in  the  process  is  to  represent  the  geometry 
with  triangular  elements  having  specific  material  characteristics.  Then  the 
distributed  mass  is  lumped  at  the  nodes  (element  corners),  and  initial  veloci- 
ties are  assigned  to  represent  the  motion  at  impact.  If  the  problem  involves 
explosive  detonation,  there  may  be  no  initial  velocities  and  the  initial  condi- 
tions are  established  by  allowing  the  detonation  process  to  begin  at  a pre- 
determined point. 


AXIS  OF  SYMMITRY 


INITIAL  CONDITIONS  INT1GRA1  ION  LOOP 


Figure  1.  Schematic  Representation  of  the  Finite 
Element  Computational  Technique 
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After  the  initial  conditions  are  established,  the  integration  loop  begins  as 
shown  in  Figure  1.  The  first  step  is  to  obtain  displacements  and  velocities 
of  the  nodes.  If  it  is  assumed  the  lines  connecting  the  nodes  (element  edges) 
remain  straight,  then  the  displacements  and  velocities  within  the  elements 
must  vary  linearly.  From  these  displacements  and  velocities,  the  strains 
and  strain  rates  within  the  elements  can  be  obtained.  Since  the  strains  and 
strain  rates  are  generally  derivatives  of  linear  displacement  and  velocity 
functions  within  the  elements,  the  resulting  strains  and  strain  rates  are 
generally  constant  within  the  elements. 

The  stresses  in  the  elements  are  determined  from  the  strains,  strain  rates, 
internal  energies  and  material  properties.  Since  the  strains  and  strain 
rates  are  constant  within  the  elements,  the  stresses  are  also  constant.  The 
stresses  are  obtained  by  combining  elastic  or  plastic  deviator  stresses  with 
hydrostatic  pressure,  Hie  deviator  stresses  represent  the  shear-strength 
capability  of  the  material,  and  the  hydrostatic  pressure  is  obtained  from  the 
volumetric  strain  and  internal  energy  of  the  element.  When  the  hydrostatic 
pressure  is  obtained  for  an  explosive  element,  the  pressure  can  be  activated 
at  a predetermined  time  (based  on  the  detonation  velocity)  or  when  the  element 
is  sufficiently  compressed.  An  artificial  viscosity  is  also  included  to  damp 
out  localized  oscillations  caused  by  representing  continuous  media  with 
lumped  masses. 

After  the  element  stresses  are  determined,  it  is  necessary  to  obtain  concen- 
trated forces  at  the  nodes.  These  forces  are  statically  equivalent  to  the  dis- 
tributed stresses  within  the  element  and  are  dependent  on  the  element  geom- 
etry and  the  magnitude  of  the  stresses.  When  the  concentrated  forces  are 
applied  to  the  concentrated  masses,  the  nodal  accelerations  are  defined,  and 
the  equations  of  motion  are  applied  to  determine  new  displacements  and 
velocities.  The  integration  loop  is  then  repeated  until  the  time  of  interest 
has  elapsed. 
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Another  feature  of  the  basic  technique  is  the  ability  to  represent  sliding  be- 
tween two  surfaces.  This  is  accomplished  with  a momentum  exchange 
principle  which  allows  for  closing,  sliding  and  separation  of  two  surfaces. 

It  should  be  noted  that  the  integration  time  increment  must  be  properly  con- 
trolled to  prevent  numerical  instability.  This  is  accomplished  by  limiting 
the  time  increment  to  a fraction  of  the  time  required  to  travel  across  the 
minimum  altitude  of  the  element  at  the  sound  velocity  of  the  material.  This 
also  ensures  that  the  time  increment  is  less  than  the  shortest  period  of 
vibration  of  the  system. 


2.  1 GEOMETRY 


A typical  axi-symmetric  triangular  element  is  shown  in  Figure  2.  It  is  geo- 
metrically defined  by  nodes  i,  j and  m,  with  the  mass  of  the  element  being 
equally  distributed  to  concentrated  masses  at  the  nodes.  When  a node  is 
contained  by  more  than  one  element,  the  total  mass  at  node  i,  M^,  is  equal 
to  1/3  the  mass  of  all  elements  which  contain  that  node.  For  an  axi-symmetr 
finite  element  model,  the  concentrated  masses  can  be  visualized  as  concen- 
tric circular  rings  contained  in  planes  which  are  perpendicular  to  the  axis  of 
revolution.  These  rings  can  move  up  and  down  along  the  axial  direction 
(Z  axis)  and  they  can  expand  and  contract  in  the  radial  direction  (R  axis). 

In  addition,  they  are  allowed  to  experience  rotations  Q about  the  axis  of 
revolution.  The  coordinates  of  node  i are  designated  r.,  z.,  (F,  and  the 
radial,  axial  and  tangential  velocities  are  designated  ir,  v.,  w.,  where 


2.2  STRAINS  AND  STRAIN  RATES 

"Die  incremental  strains  which  occur  during  each  cycle  of  integration  are  ob- 
tained by  multiplying  the  strain  rates  by  the  integration  time  increment. 

The  strain  rates  are  obtained  from  the  current  geometry  of  the  element  and 
the  velocities  of  the  nodes.  If  it  is  assumed  that  the  lines  connecting  the 
nodes  remain  straight,  the  displacements  and  velocities  within  each  element 
must  vary  in  a linear  manner.  Then  the  velocities  within  the  element  can  be 
expressed  as 

ti  = a1  + a2  r + a^z  (1) 
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v - a4  + o5r  + a6 z 


(2) 
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i 
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« 
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w = a7  + a8r  + a0z  (3) 

where  . . . ag  are  geometry  and  velocity-dependent  constants.  It  is  pos- 
sible to  solve  for  Oj,  a^,  by  substituting  the  radial  velocities  and  coordi- 
nates of  nodes  i,  j,  m into  Equation  (1).  This  gives  three  equations  and 
three  unknowns  such  that  Equation  (1)  can  be  expressed  in  terms  of  the  ele- 
ment geometry  and  nodal  velocities. 

“ = SC[(,i+bir'K:iI>ii  + (aj+bjr+CjZ,ij  + <am+bmr+cmI,ilm]  (4) 

where  a.  = r.z^  - r z.f  b.  = z.  - z . c.  = r - r.,  and  A is  the  cross- 
ljmmjijmimj 

sectional  area  of  the  element  in  the  R-Z  plane.  The  other  velocities  (v,  w) 
are  identical  to  Equation  (4)  except  the  radial  velocities  are  replaced  by  the 
axial  and  tangential  velocities. 

After  the  velocities  are  obtained,  it  is  possible  to  determine  the  normal 

strain  *ates  (e^,  e^,  e0),  the  shear  strain  rates  (Yrz,  Ytq.  yz0)  and  the 

localized  rotational  spin  rate  of  the  element  in  the  R-Z  plane  (u  ). 

rz 


• 3u 

Cr  = 3? 

(5) 

e = -2^ 
ez  3z 

(6) 

%-4 

(7) 

y = + 

'rz  dz  3r 

(8) 

y = - E 

rrQ  or  r 

(9) 
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i 

i j 
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z 0 


u_ 


rz 


dw 

ST 

(10) 

1 /dv  Bu'N 

2"Vdr  dz ) 

(11) 

It  can  be  seen  that  Equations  (5),  (6),  (8),  (10)  and  (11)  are  derivatives  of 
linear  functions  and  are  therefore  constant  within  the  element.  Equations 
(7)  and  (9)  involve  averages  of  the  nodal  velocities  (u  and  w)  and  the  three 
radii  (r),  so  they  are  not  necessarily  constant.  It  is  also  necessary  to  use 
an  equivalent  strain  rate  which  is  expressed  as 


i ■ Vt  [(  i - ;z»2  + <*r  - + + T frrz2 


(12) 


+ y 2 + v 2)] 

7r0  + yz6  ' J 


An  equivalent  plastic  strain,  e , is  then  obtained  by  integrating  e with  re- 
spect to  time  during  plastic  flow. 


- t+At 

C 

P 


e * + eAt 
P 


(13) 


where  At  is  the  integration  time  increment.  It  should  also  be  noted  that  sub 

sequent  computations  will  involve  deviator  strain  rates  (e  , e , e„)  which 

. . r . z 0 

are  readily  obtained  from  the  normal  strain  rates  (e^,  e^). 


2.  3 STRESSES  AND  PRESSURES 

The  stresses  in  the  elements  are  determined  from  the  strains,  strain  rates, 

internal  energies  and  material  properties.  The  three  normal  stresses  (o  , 

r 

a , afl)  are  expressed  in  terms  of  deviator  stresses  (s  , s , s ) hydrostatic 
pressure,  P,  and  artifical  viscosity  Q. 
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°r  * sr  “ {P  + Q> 


°z  c 8z  - (P  + 


ae  = Bo  " (p  + Q) 


Trial  values  of  the  deviator  stresses  at  time  * t + At  are 


t+At 
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= s1 

+ 

2Ge 

At 

- 2t*  u>  At 
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r 

r 

rz  rz 

t+At 

i 

= s* 
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to 

At 

+ 2Tt  u At 
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rz  rz 

t+At 
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II 

CD 

+ 

2G®0 

At 

In  Equation  (17)  the  first  term  (s  ) is  the  radial  stress  at  the  previous  time 

• 1 

and  the  second  term  (2GerAt)  is  the  incremental  stress  due  to  the  incremental 
strain  (e  At)  during  that  time  increment,  where  G is  the  elastic  shear  mod- 

r t 

ulus.  The  third  term  (2t  u At)  is  due  to  shear  stresses  from  the  previous 

rz  rz 

time  increment,  which  now  act  as  normal  stresses  due  to  the  new  orientation 
of  the  element  caused  by  an  incremental  rotation  («rzAt)  during  the  time  in- 
crement. The  axial  stress  has  the  same  form  as  the  radial  stress,  and  tan- 
gential stress  is  also  similar  except  there  is  no  contribution  from  rotated 
shear  stresses. 


The  trial  values  of  the  shear  stresses  are  formulated  in  similar  manner. 

' Trz  + GVzAt  + <ar  ' <4>  “rz"  <20> 


■ Tr0  + 


t+At  t _ • 
rz  0 = Tz  0 + Gyz6At 


The  rotational  terms  can  only  be  significant  for  stresses  in  the  R-Z  plane 
since  the  elements  can  rotate  without  experiencing  strain  rates.  For  the 
other  stresses  involving  9,  the  analogous  rotations  cannot  occur  without  in- 
ducing strain  rates,  and  the  effect  of  the  rotational-induced  stresses  is  small 
compared  to  the  strain  rate  induced  stresses. 

It  should  ]je  recalled  that  Equations  (17)  through  (22)  represent  trial  values 
of  the  stresses  and  they  may  need  to  be  reduced  if  they  violate  the  Von  Mises 
yield  criterion.  An  equivalent  stress  is  given  by 

’ " Vl<Sr  + al  + 80>  + 3(Trz  + Tr0  + ri0>  <23) 

If  cf  is  not  greater  than  the  equivalent  tensile  strength  of  the  material,  S, 
the  final  deviator  and  shear  stresses  are  as  given  in  Equations  (17)  through 
\22).  If  5 is  greater  than  5,  then  the  stresses  in  Equations  (17)  through  (22) 
should  be  multiplied  by  the  factor  (3/tf).  When  the  reduced  deviator  and  shear 
stresses  are  put  into  Equation  (23),  the  result  is  always  = 3. 

The  equivalent  tensile  strength  of  the  material  may  be  dependent  on  many 
factors,  including  strain,  strain  rate,  pressure  and  temperature.  It  is 
well-known  that  many  materials  behave  differently  under  dynamic  impact 
than  under  static  testing  conditions.  With  few  exceptions,  however,  a pre- 
cise definition  of  material  behavior  under  dynamic  conditions  is  not  available. 
There  is  also  much  to  be  learned  about  fracture  characteristics  under  these 
dynamic  conditions.  The  current  version  of  EPIC -2  allows  the  equivalent 
tensile  strength  to  be  determined  from 

S = Sg  [1  +Cj  log  (e)  ] [1  +C2P  +C3P2]  [C4+Cf)T]  (24) 

In  this  formulation,  Sg^  is  generally  taken  to  be  the  static  stress,  which  is 
dependent  on  the  equivalent  plastic  strain  of  Equation  (13).  The  three  brack- 
eted terms  allow  the  static  stress  to  be  altered,  based  on  strain  rate, 


10 


F 
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pressure  and  temperature.  If  the  constants  are  Cj  =C2  BCj  = C&  = 0 and 
C4  *1.0,  then  I?  * S f , which  is  the  strain-dependent  static  stress.  The  first 
bracketed  term  can  increase  the  strength  due  to  the  equivalent  strain  rate, 
e,  described  in  Equation  (12).  In  the  second  bracketed  term,  P * P + Q,  so 
the  effect  of  this  term  is  to  increase  the  strength  due  to  pressure.  The  final 
bracketed  term  includes  the  temperature,  T,  of  the  element,  which  is  ob- 
tained from 


E 

T = T (25) 

O C p 

sro 

where  Tq  is  the  initial  temperature  of  the  element,  Eg  is  the  internal  energy 

per  original  unit  volume  (defined  later),  C is  the  specific  heat  and  p is  the 

s o 

initial  density. 

Material  fracture  is  currently  dependent  on  the  equivalent  plastic  strain  7 , 
(Equation  13)  and  the  volumetric  strain,  e = V/Vo-l,  where  V and  Vo  re- 
present the  current  and  initial  volumes  of  the  element.  When  the  fracture 
criterion  has  been  met,  the  equivalent  tensile  stress  is  set  equal  to  zero, 
so  no  shear  stresses  can  be  developed  in  the  failed  element.  Likewise,  no 
tensile  stresses  are  allowed  to  develop.  The  net  result  is  that  a failed  ele- 
ment tends  to  act  like  a liquid  inasmuch  as  it  can  develop  only  hydrostatic 
compression  with  no  shear  or  tensile  stresses.  Another  option  is  also 
available  which  sets  all  stresses  (including  the  pressure)  equal  to  zero. 

The  hydrostatic  pressure  is  dependent  on  the  volumetric  strain  and  the  in- 
ternal energy  in  the  element.  The  EPIC -2  code  uses  the  Mie-Grilneisen 
equation  of  state  in  the  general  form  P = Pv  + rEg  (1+p),  with  the  complete 
expression  given  by 

P = (KjP  + K2m2  + K3m3)  (1  --^-)  + TEg  (1  +p)  (26) 
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where  p = Vo/V-1,  Kj,  K and  Kg  are  material-dependent  constants  and  r 
is  the  GrUneisen  coefficient. 

Since  the  pressure  can  be  significantly  affected  by  the  internal  energy  E , it 

S 

is  desirable  to  solve  the  pressure  and  energy  equations  simultaneously.  This 
gives 


t+At  E1  - • 5 L(P  + Q)*  + Qt+At  + Pt+At]  e At  + AE 

E = _ v v (27) 

S 1 + .5T(1  +p)  evAt 

where  is  the  volumetric  strain  rate  and  AE^  is  the  internal  energy  gen 
erated  by  the  deviator  and  shear  stresses  during  the  previous  cycle. 


AE  , = (S  e + IT  e + Enen  + T y + f .y  . + f ,y  .) 
d r r z z GO  rz'rz  rG'rG  zG'zO 


(V/V  ) At 
o 


(28) 


The  "bars"  on  the  deviator  and  shear  stresses,  and  the  volume,  represent 
averages  of  these  values  at  times  t and  t+At.  After  the  total  internal  energy 
at  time  t+At  has  been  determined  from  Equation  (27),  the  pressure  at  time 
t+At  is  determined  from  Equation  (2G). 


For  explosive  detonation,  the  hydrostatic  pressure  is  determined  by  a pro- 
cedure similar  to  that  used  in  the  HEMP  code.*  The  pressure,  determined 
by  the  Gamma  Law,  is  given  by 

P = F (y  - 1)  Eg  (1  +p)  (29) 

where  F is  the  burn  fraction  (0  sF  s 1,0),  y is  a material  constant,  and  E 

s 

is  the  internal  energy  per  initial  unit  volume.  Unlike  the  solid  and  liquid 
materials,  the  explosive  contains  internal  energy  in  the  initial  condition. 
The  pressure  and  energy  equations  are  solved  simultaneously  in  a manner 


. jy-sawf1 


similar  to  that  of  Equations  (26)  and  (27).  The  explosive  is  effectively  ini- 
tiated with  the  burn  fraction,  which  is  dependent  on  the  time  for  detonation 
to  arrive  and  travel  through  the  element,  or  the  compressed  state  of  the  ele 
merit.  The  burn  fractions  for  these  two  conditions  are 


F 


F 


(t8  - t)  D 

Wo 

1 - V/Vo 


(30a) 

(30b) 


In  Equation  30a,  A is  the  initial  area  of  the  element  and  t is  the  time  re- 
o s 

quired  for  the  detonation  wave  to  reach  the  node  closest  to  the  point  of  ini- 
tiation when  traveling  at  the  detonation  velocity,  D.  This  formulation  tends 
to  spread  the  wave  front  over  a limited  number  of  elements.  Equation  30b 
gives  the  burn  fraction  in  terms  of  the  compressed  state,  where  V,-..  = 

V ^*>0 

+ ^ is  the  Chapman-Jouquet  relative  volume.  This  allows  a converging 
detonation  wave  to  travel  at  a velocity  greater  than  D.  The  minimum  and 
maximum  allowable  values  of  F are  0.  and  1.  0. 


The  artificial  viscosity  is  combined  with  the  normal  stresses  to  damp  out 
localized  oscillations  of  the  concentrated  masses.  It  tends  to  eliminate  spur- 
ious oscillations  which  would  otherwise  occur  for  wave  propagation  problems. 

2 

This  technique  was  originally  proposed  by  Von  Neumann  and  Richtmyer  and 

3 

has  been  expanded  for  use  in  various  computer  codes.  It  is  expressed  in 
terms  of  linear  and  quadratic  components  and  is  applied  only  when  the 
volumetric  strain  rate  is  negative: 

Q = CLpcsh  Uvl  + c2  Ph2(ev)2  for  < 0 
Q = 0 for  ev  2 0 


(31) 


TP"’ 


E I 


where  c and  p are  the  sound  velocity  and  density  of  the  material  and  h is  the 
s 

minimum  altitude  of  the  triangle.  Typical  values  used  for  the  dimensionless 

2 1 

coefficients  are  C T =0.5  and  C =4.0. 

Li  O 


The  sound  velocity  for  solid  and  liquid  elements*  is  obtained  from 

c2  = -1—  [K  (1  - Til)  + K.(2|I  - 1.5rw2)  + K„  (3u2  - 2rW3) 
s Po  1 2 A 

+ r es  + rP/U  + p)] 

and  the  sound  velocity  for  explosive  elements  is  obtained  from 


yP 


(32) 


(33) 


The  sound  velocities  are  also  used  for  determination  of  the  integration  time 
increment. 


2.4  CONCENTRATED  FORCES 

After  the  element  stresses  are  obtained,  it  is  necessary  to  determine  con- 
centrated forces  to  act  on  the  concentrated  mass  at  the  nodes.  This  is  done 
by  obtaining  the  concentrated  forces  which  are  statically  equivalent  to  the 
distributed  stresses  in  the  elements.  The  radial,  axial  and  tangential  forces 
acting  on  node  i of  an  element,  are 


F*  = -ttT  [(z.  - z )cf  + (r  - t.)t  l-^-rrAa 
r L j m r m j rz  J 3 6 


F1  = -tt?  T(r  - r.)a„  + (z.  - z )t  1 
z m j z ] m rzJ 


(34) 

(35) 


14 


1 


F'  = -iff  Cf  (*j  - zm)rr6  + (rm  - rj)Tie] 


In  Equation  (36)  the  factor  T/rj  is  included  in  the  expression  Involving  the 
shear  stress  in  the  R-0  plane.  This  is  necessary  to  satisfy  the  equilibrium 
condition  in  this  plane  which  is  -f  2Tffl  = 0.  5 It  is  clear  from  this  con- 
dition that  the  shear  stress  cannot  be  a non-zero  constant,  but  rather  must 

vary  as  a function  of  the  radius.  This  stress  is  obtained  from  y in  Equa- 

rt/ 

tion  (9)  and  is  not  necessarily  constant  since  it  involves  averages  of  the  three 

nodal  displacements  and  radii.  Therefore,  the  factor  T'/r^  in  Equation  (36) 

represents  the  effect  of  a variable-shear  stress  in  the  element.  The  net 

forces  at  node  i (F* , f\  f!.)  are  the  sum  of  the  forces  from  each  of  the  in- 
r t o 

dividual  element  forces  at  that  node. 


2.  5 EQUATIONS  OF  MOTION 

The  equations  of  motion  are  integrated  by  assuming  a constant  velocity  for 
each  time  increment.  The  acceleration  of  node  i in  the  radial  direction  at 
time  = t is 


= — L+r*  (et_)2 
1 V V / 


The  first  term  is  due  to  stress-induced  forces  and  the  second  term  is  due  to 
spin.  The  updated  velocity  for  the  next  time  increment  is 


u*+  = (u*~  + u*  5f)  (1  - At/r|) 


where  uj  is  the  constant  velocity  for  the  previous  time  increment  and  51  is 
the  average  of  the  two  integration  time  increments  about  time  = t.  The  ex- 
pression in  the  second  set  of  parentheses  can  be  used  to  damp  out  the  radial 


velocities  to  give  steady-state  solutions  for  spinning  bodies.  If  the  constant, 
Cjy  is  set  equal  to  twice  the  sound  velocity  of  the  material,  the  system  will 
be  approximately  critically  damped  and  the  steady-state  solution  will  be 
rapidly  attained.  This  will  later  be  demonstrated  in  an  example.  Finally, 
the  new  radial  displacement  at  time  = t+At  is 

u|+At  = u*  + u*+  At  (39) 

The  equations  of  motion  for  the  axial  direction  have  a similar  form  except 
the  spin  effects  for  the  acceleration  and  the  damping  effects  for  the  velocity 
are  not  included. 


For  the  circumferential  equations  of  motion  it  is  necessary  to  consider  the 
angular  momentum,  H,  of  each  concentrated  mass.  This  gives 


Hi+  = Hi"  + r\veM  (40) 

By  substituting  H.  = 0 . r.  into  Equation  (40),  it  is  possible  to  determine 
the  updated  rotation  velocity. 


r]  5F 

rETZE? 

(ri  } 


(41) 


It  should  be  noted  that  even  if  the  net  circumferential  force,  is  equal  to 
zero,  it  is  possible  for  the  spin  to  change  if  the  radius  changes  between  times 
t and  t+At.  It  is  therefore  necessary  to  obtain  the  new  radial  position  at 
t+At  before  obtaining  the  new  rotation  velocity  at  t+At. 


The  integration  time  increment  must  be  controlled  to  prevent  numerical  in- 
stability. This  is  done  by  limiting  the  time  increment  to 


At  - C, 


2 2 

where  g = Cq  Q/p,  hmjn  is  the  minimum  altitude  of  the  triangle  and  cg  is 
the  sound  velocity.  ^ The  constant,  C^,  must  be  less  than  unity  to  ensure 
that  At  is  always  less  than  the  time  required  to  travel  across  the  shortest 
dimension  of  the  triangle  at  the  sound  velocity  of  the  material.  This  criterion 

g 

also  ensures  that  At  is  leps  than  the  lowest  period  of  vibration  of  the  system. 
The  EPIC -2  program  restricts  the  time  increment  from  increasing  more 
than  10  percent  per  cycle. 


2.  6 SUDINC.  SURFACES 

It  is  sometimes  necessary  to  allow  for  sliding  to  occur  between  two  surfaces. 
The  important  steps  for  the  sliding-surface  technique  are  summarized  as 
follows: 

• Identify  a "master"  sliding  surface  defined  by  a specified  row 
of  master  nodes. 

• Identify  a "slave"  surface  (or  region)  defined  by  a specific  row 
(or  group)  of  slave  nodes.  The  slave  node  spacing  should  not 
be  significantly  greater  than  the  master  node  spacing  since 
this  could  introduce  localized  deformations  in  the  master  sur- 
face at  the  slave  node  locations. 

• For  each  integration  time  increment,  apply  the  equations  of 
motion  to  both  the  master  nodes  and  the  slave  nodes  in  the 
usual  manner. 


For  each  slave  node,  find  the  portion  of  the  master  surface 
(defined  by  two  master  nodes)  whose  projection  contains  that 
slave  node. 
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• Check  to  determine  if  there  is  interference  between  the  slave 
node  and  the  master  surface. 

• If  there  is  interference,  place  the  slave  node  on  the  master 
surface  in  a direction  normal  to  the  master  surface. 

• Determine  the  translational  momentum  of  the  two  master 
nodes  and  the  slave  node  in  a direction  normal  to  the  master 
surface.  Also  determine  the  angular  momentum  of  these 
three  nodes  about  a point  on  a line  which  contains  the  two 
master  nodes. 

• Determine  updated  velocities  of  the  two  master  nodes  and  the 
slave  node  in  a direction  normal  to  the  master  surface. 

• Adjust  the  volumes  (optional)  of  the  slave  elements  which  have 
a master  node  located  between  the  two  slave  nodes  of  the  slave 
element. 

This  technique  is  illustrated  in  Figure  3.  It  can  be  seen  that  slave  node  s 
travels  from  A (at  time  = t - At)  to  B (at  time  = t)  using  the  equations  of  mo- 
tion from  subsection  2.  5.  Since  the  velocity  is  constant  for  each  time  in- 
crement, the  node  travels  in  a straight  line  as  shown.  The  master  surface 
is  also  shown  for  time  = t,  and  it  can  be  seen  that  slave  node  s,  at  position 
B,  has  passed  through  the  line  connecting  master  nodes  i and  j.  The  slave 
node  is  next  moved  to  position  C,  which  is  on  the  master  surface.  The  line 
from  B to  C is  normal  to  the  master  surface  and  therefore  affects  the  equa- 
tions of  motion  only  in  this  normal  direction.  This  also  has  the  effect  of  as- 
suming a frictionless  surface.  ^ 

• 
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MASTER  MATERIAL  AT  TIME  ■ t 


Figure  3.  Sliding  Surface  Procedure 


The  next  step  consists  of  updating  the  normal  velocities  (V..  V.,  Vg)  of  the 
three  nodes.  Since  three  unknowns  are  involved,  three  conditions  must  be 
specified  to  obtain  a solution.  Two  of  these  conditions  are  that  the  transla- 
tional and  rotational  momenta  of  the  three  nodes  should  be  conserved.  The 
third  condition  is  that  the  normal  velocity  of  the  slave  node  is  equal  to  the 
normal  velocity  of  the  master  surface  at  the  location  of  the  slave  node.  The 
resulting  three  equations  are 


W.V.t+r.'  + M.Vt+r.'  + M V t+r  ' = M.V.t_r'  + M.V^V 


+ M V t_r' 
s s s 


vt+  + (Vt+  - Vt+) 
i 3 i 


_/  „ / 
r - r . 
s 1 

rf  - r! 
3 i 


where  the  three  unknowns  are  the  updated  velocities  (V*+,  vV*\  V*+).  The 

1 ] s 

coordinates  in  Equations  (44)  and  (45)  (r! , r r *)  are  along  the  R ' axis 

* J ® 

which  includes  the  two  master  nodes.  It  should  be  noted  that  the  use  of  Equa- 
tion (45)  makes  the  solution  slightly  dependent  on  the  order  that  the  slave 

7 

nodes  are  processed.  A previously  used  procedure  was  to  transfer  the  mo- 
mentum change  imposed  on  the  slave  node  (determined  from  the  movement 
from  B to  C in  Figure  3)  to  the  adjacent  master  nodes.  The  use  of  Equation 
(45)  was  selected  to  reduce  the  "rattling"  on  the  surface  by  better  matching 
the  normal  velocities.  An  exact  velocity  match  at  all  slave  node  locations 
could  be  obtained  by  repeatedly  applying  Equations  (43)  to  (45)  until  the  V*+ 
velocities  were  equal  to  the  V*  velocities.  This  would  eliminate  any  order 
dependency.  It  should  be  emphasized  that  the  order  dependency  is  very  min- 
imal since  the  net  momenta  are  conserved  for  the  two  master  nodes  and  the 
slave  node.  It  should  also  be  noted  that  small  errors  are  introduced  into  the 
center  of  gravity  (eg)  positions  when  the  slave  node  is  moved  to  the  master 
surface,  since  there  is  no  corresponding  movement  of  the  master  nodes, 
which  receive  only  instantaneous  velocity  changes.  The  eg  errors  are  very 
small,  however,  and  do  not  significantly  affect  the  solutions. 


2.  7 SEVERE  DISTORTIONS 

A characteristic  of  a triangular  formulation  is  that  it  is  better  suited  to  re- 
present severe  distortions  than  is  a quadrilateral  formulation.  This  is  due 
to  a triangle's  resistance  to  allowing  a node  to  cross  the  opposite  side  of  the 


triangle  during  the  high  distortion.  It  is  apparent  that  the  cross- sectional 
area  of  a triangle  approaches  zero  as  a node  approaches  the  opposite  side  of 
that  triangle.  Since  the  hydrostatic  pressure  of  Equations  (26)  and  (29)  is 
inversely  proportional  to  the  volume  of  the  element,  as  the  volume  approaches 
zero  the  pressure  approaches  infinity.  As  a result,  there  is  a very  large 
resistance  to  this  mode  of  distortion  and  the  triangle  will  not  go  through  a 
zero  volume  configuration  provided  the  time  increment  is  properly  controlled. 

The  quadrilateral  formulation  does  not  have  this  desirable  characteristic. 
Figure  4 shows  three  quadrilateral  configurations  which  have  equal  cross- 
sectional  areas,  Aq.  Configuration  A shows  a rectangular  cross-section 
with  nodes  1,  j,  m and  p defining  the  initial  geometry'.  Configuration  B re- 
presents a deformed  element  with  nodes  j and  m occupying  the  same  posi- 
tion. A triangular  element  with  two  nodes  occupying  the  same  position 
would  result  in  a zero  cross-sectional  area.  For  the  quadrilateral  element, 
however,  nodes  i and  p can  be  displaced  to  allow  the  initial  area  to  be  re- 
tained. Configuration  C shows  overlapping  nodes  and  the  net  area  of  the 
quadrilateral  is  equal  to  the  difference  of  the  two  triangular  cross  sectiona 
Even  this  configuration  can  exist  with  the  initial  area.  When  this  type  of 
distortion  is  achieved,  however,  the  formulation  generally  breaks  down  and 
numerical  instability  often  occurs.  The  deformations  represented  by  Con- 
figurations B and  C usually  occur  when  the  grid  is  highly  distorted.  These 
deformations  are  resisted  by  the  deviator  stresses,  but  the  deviator  stresses 
are  limited  by  the  strength  of  the  material  and  do  not  provide  the  potential 
magnitude  of  resistance  as  does  the  hydrostatic  pressure  for  the  triangular 
elements. 
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CONFIGURATION  A 


Figure  4.  Possible  Distortions  of  a Quadrilateral  Element 


2.  8 BASIC  VERIFICATION  EXAMPLES 


Several  basic  verification  examples  are  included  to  illustrate  and  verify  the 
computational  technique.  The  first  examples  are  shown  in  Figure  5 and  they 
demonstrate  the  capability  of  the  technique  to  provide  accurate  compressive 
wave  propagation  solutions.  Figure  5a  shows  a radially  restrained  column 
of  water  impacting  a rigid  wall  at  300m /s  (Mach  0.  2).  Even  though  this  re- 
duces to  a one-dimensional  problem  it  is  a severe  test  since  the  water  is 
significantly  compressed  and  the  shock  velocity  is  a function  of  the  impact 
velocity,  and  therefore  higher  than  the  acoustic  velocity.  The  dimensionless 
pressure  is  defined  as  P*  = P/p  C V , where  P,  p , C and  V represent 
the  pressure,  initial  density,  acoustic  velocity  and  impact  velocity,  respec- 
tively. The  pressure  is  shown  at  a time  shortly  before  the  wave  reaches  the 
far  end  of  the  column.  It  can  be  seen  that  the  finite  element  solution  with 
artificial  viscosity  agrees  with  the  analytic  solution  along  the  length  of  the 
column,  with  only  slight  differences  occurring  at  the  leading  edge  of  the  wave. 
The  solution  obtained  without  artificial  viscosity  is  also  shown  and  it  can  be 
seen  that  the  pressure  radically  oscillates  about  the  analytic  solution. 

A similar  example  is  shown  in  Figure  5b  where  a wave  propagation  solution 
iu  demonstrated  for  a detonated  explosive.  Here  the  dimensionless  pressure 
is  defined  as  P'v  3 P/[p^D“/(y  +1)  ] where  P,  p and  D are  the  pressure,  ini- 
tial density  ami  detonation  velocity,  respectively,  and  y is  a material-depen- 
dent constant.  Again  it  can  be  seen  that  there  is  excellent  agreement  between 
the  finite  element  and  analytic  solutions. 

Another  compressive  wave  propagation  example  is  shown  in  Figure  6.  This 
example  shows  the  effect  of  material  strength  with  a higher-velocity  elastic 
wave  at  the  leading  edge  of  the  shock  front  and  the  elastic  unloading  at  the 
rear.  This  solution  is  in  good  general  agreement  with  that  presented  in  Re- 
ference 1. 
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COLUMN  Of  HATER  WITH  RADIAL  Rf STRAIN! 
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KICID  ► / 1 INITIAL  GEOMEU 

boundary  .luvninuuuuuuninnnnnunnnmvvuJ  at  impact  or 

i * V„  • JOOM/S 

1 ♦ FINITl  ELEMENT  MOOLI  0 

3?ctui mnniui mumwvwwi mini mtu id  DI5PlACcl“t 

AT  3-  • O.SB 

*r  i 

I IMPACT  EXAMPLE  J 

LV  jyj  | 

lj}f  V | | |!  M — ANALYTIC  SOLUTION 

S ) FINITE  lltMtNI  S01U1I0N 

3 n II  ! 1 UUTU  ADTinriAl  uicrnci 


analytic  SOLUTION 


* h i V b ' i.  i 1 1 j I 


FINITE  ULMENI  SOLUTION 
W ITH  artificial  viscosity 


ao-o-oo  F INITL  ELEMENT  SOLUTION 

WITHOUT  ARTIFICIAL  VISCOSITY 


I I !! 


:iW 


\ 

»♦*♦♦♦ 


DETONATION  I 
INITIAILD  AT  *■  } 
RIO  10  BOUNDARY  ! 


LLHUMN  Of  EXPLOSIVE  WITH  RADIAL  RESTRAINT 

IN  * “T  INITIAl  GtOMtlRI 

at  ► u umuin : iinimniuuiuiiiniuii i i ; i : i lixJ  at  ouonahon 

3ARY  ! 

.li  1 1 1 1 1 ; n m 1 1 1 1 u : 1 1 1 ii  in  unir.mmniun  il  “Ils,’oit'011cJ,1 
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Figure  5.  Wave  Propagation  Due  to  Impact  and  Explosive  Detonation 


ran 


| IMPACT  VUOCITY  • 800  M/S 


BAR  1 j BAR  2 


Figure  6.  Wave  Propagation  Resulting  From  the  Impact  of  two  Bars 


The  capability  to  provide  shear  wave  propagation  solutions  is  shown  in  Figure 
7 where  a torque  is  suddenly  applied  to  the  end  of  a thin-walled  cylinder.  It 
can  be  seen  that  there  is  good  general  agreement  between  the  two  solutions, 
although  the  EPIC -2  solution  does  oscillate  about  the  analytic  solution.  These 
oscillations  exist  because  there  is  no  artificial  viscosity  for  the  shear  mode 
of  deformation.  The  previously  defined  artificial  viscosity  in  Equation  (31) 
is  affected  only  by  volumetric  strain  rates.  The  need  for  artificial  viscosity 
in  the  shear  mode  is  not  as  critical  as  for  the  compressive  mode,  however. 
This  is  because  the  magnitude  of  the  shear  stress  is  boui  led  by  the  strength 
of  the  material  whereas  there  is  no  upper  bound  for  the  magnitude  of  com- 
pressive waves.  With  the  unbounded  compressive  stress  it  can  be  seen  in 
Figure  5a  that  very  radical  oscillations  can  occur  if  the  artificial  viscosity 
is  omitted.  Since  these  radial  oscillations  cannot  occur  in  the  shear  mode, 
due  to  the  limitation  on  the  strength  of  the  material,  the  accuracy  of  the 


25 


APPLIED  TORQUE 


c 


THIN  WALLED  CYLINDER  (4*  10) 


R 


± 

T 


H 


Figure  7.  Shear  Wave  Propagation  in  a Hollow  Cylinder 

solution  in  Figure  7 is  generally  adequate  and  no  artificial  viscosity  is  in- 
cluded for  the  shear  mode  of  deformation. 

A characteristic  of  spinning  bodies  is  that  transient  dynamic  responses  are 
introduced  if  the  spin  is  suddenly  applied.  This  can  be  illustrated  by  consid- 
ering the  case  of  a hoop  subjected  to  a suddenly  applied  spin.  It  can  be  shown 
that  the  radial  stiffness  of  a ring  is  K = 2 ttAE/R  and  the  mass  is  M = 2 

ttR  Ap  where  A is  the  cross-sectional  area,  E is  the  modulus  of  elasticity, 
o 

Ro  is  the  radius  and  p is  the  density.  The  period  of  vibration  then  becomes 

T = 2ttVM/K  = 2rrR  /c  where  c = vAn/p  is  the  sound  velocity.  Figure  8 
OS  s 

shows  the  response  of  a ring  subjected  to  a suddenly  applied  spin  u.  For 
this  problem  the  hoop  is  represented  by  three  nodes  and  one  triangular  ele- 
ment. It  can  be  seen  that  the  undamped  spinning  hoop  behaves  like  a single- 
degree -of-freedom  oscillator  subjected  to  a suddenly  applied  force.  The  re- 
sponse is  plotted  as  dimensionless  stress  as  a function  of  dimensionless  time. 
For  the  case  of  the  spinning  hoop  the  dimensionless  stress  is  the  ratio  of  the 
dynamic  stress  to  the  steady-state  stress.  Likewise,  the  dimensionless 
time  is  the  ratio  of  elapsed  time  to  the  period  of  vibration  of  the  spinning 


26 


hoop.  Both  the  magnitude  of  the  response  and  the  period  of  vibration  are 
identical  to  those  of  a s ingle-degree -of-freedom  oscillator. 


A 


SPINNING  HOOP 
- UNDAMPED 


SPINNING  DISK 


UNDAMPED 


DIMENSIONLESS  TIME  (j^) 


Figure  8.  Dynamic  Response  of  a Spinning  Hoop  and  a Spinning  Disk 
Subjected  to  a Suddenly  Applied  Spin 


Since  a spinning  projectile  probably  attains  a steady-state  condition  prior 
to  striking  the  target,  it  is  desirable  to  represent  this  steady-state  condi- 
tion with  the  computer  simulation.  This  can  be  done  by  assuming  the  system 
is  critically  damped  until  a steady-state  condition  is  achieved.  Referring 
again  to  a single-degree-of-free  system,  the  critical  damping  coefficient  is 
CCRIT  = 2Mcs/Ro-  Since  the  resulting  radial  force  is  FCRIT  = -CCRIT  R. 
the  resulting  velocity  change  during  a time  increment  is  All  = 

= -2cfl2lk/Ro.  This  provides  the  basis  for  the  damping  effect  in  Equation 
(38)  where  the  constant  is  set  equal  to  twice  the  sound  velocity  (2c„). 


The  response  of  a spinning  disk  is  also  shown  in  Figure  8.  The  basic  re- 
sponse is  similar  to  that  of  the  hoop  except  both  the  period  of  vibration  and 
the  magnitude  of  the  stress  are  much  lower.  The  state  of  stress  throughout 
the  disk  is  shown  in  Figure  9 and  it  can  be  seen  that  there  is  excellent  agree- 
ment between  the  two  solutions.  Although  the  primary  function  of  the  EPIC -2 
code  is  to  obtain  solutions  for  impact  problems,  it  also  provides  a very  ef- 
ficient means  of  computing  steady-state  solutions  for  spinning  axi-symmetric 
bodies. 

Z 

^FINITE  ELEMENT  MODEL  OF  SPINNING  DISK 


Figure  9.  Stress  Distribution  in  a Spinning  Disk 

This  completes  the  discussion  of  the  basic  verification  examples.  Input  data 
for  a more  comprehensive  example  is  provided  in  the  Appendix. 
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SECTION  3 

COMPUTER  PROGRAM  DESCRIPTION 

The  EPIC -2  computer  program  contains  the  formulation  presented  in  Section 
2.  A description  of  some  of  the  characteristics  of  the  computer  program 
is  given  in  the  following  subsections. 
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3.  1 PROGRAM  ORGANIZATION  AND  SUBROUTINES 


The  organization  of  the  EPIC-2  code  is  shown  in  the  hierarchy  chart  in  Figure 
10.  It  consists  of  the  main  program.  EPIC -2,  and  18  subroutines  for  a total 
of  3140  cards.  The  EPIC -2  code  is  core  contained  and  requires  76K  storage 
on  a Honeywell  6080  computer  for  a capacity  of  1000  nodes  and  2000  elements. 
Two  plotting  programs  are  also  available  to  provide  state  plots  (geometry, 
stress,  strain,  pressure  and  velocity  fields)  at  specific  times,  and  to  plot 
various  system  parameters  (energy,  momentum,  etc.)  as  functions  of  time. 
The  following  is  a brief  description  of  the  main  program  and  subroutines 
(listed  alphabetically)  shown  in  Figure  10. 


EPIC -2  This  is  the  main  program  which  calls  subroutines  HDATA, 
MATL,  GEOM,  START  and  LOOP  (77  cards) 

BREAK  This  subroutine  checks  for  fracture  (21  cards) 

ELEG  This  subroutine  generates  a series  of  triangular  elements  (61 

cards) 

ESHAPE  This  subroutine  generates  elements  for  special  shapes  includ- 
ing rods  with  various  nose  shapes,  spheres  and  plates  (236 
cards) 
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EPIC-? 

MIN 

CALLING 

PROGRAM 


GEOM 

HDATA 
HE  BURN 

LOAD 

LOOP 


MATL 

NODEG 

NSHAPE 

RECALL 

SAVE 

SDATA 

SLIDE  1 


This  subroutine  reads,  generates  and  prints  the  initial  geo- 
metry (438  cards) 

This  subroutine  prints  all  data  as  it  is  input  (29  cards) 

This  subroutine  computes  pressures  and  energies  for  explosive 
elements  (58  cards) 

This  subroutine  computes  internal  axial  and  torsional  loads  in 
an  axi-symmetric  rod  (41  cards) 

This  subroutine  computes  the  equations  of  motion,  strains, 
strain  rates,  stresses,  pressures  and  concentrated  nodal 
forces.  It  also  prints  selected  node  and  element  data  (498 
cards) 

This  subroutine  reads  and  prints  material  data  (189  cards) 

This  subroutine  generates  a row  of  nodes  (64  cards) 

This  subroutine  generates  nodes  for  special  shapes  including 
rods  with  various  nose  shapes,  spheres  and  plates  (220  cards) 

This  subroutine  reads  data  from  Tape  2 for  a restart  run  (93 
cards) 

This  subroutine  writes  data  on  Tape  2 for  a restart  run  or  a 
state  plot  (89  cards) 

This  subroutine  computes  system  data  such  as  energy,  momen- 
tum, etc. , for  the  projectile,  target  and  total  system.  It  also 
writes  this  data  on  Tape  3 for  time  plots  (157  cards) 

This  subroutine  computes  the  equations  of  motion  for  the  slave 
nodes.  If  there  is  interference  with  the  master  surface,  it 
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places  the  slave  node  on  the  master  surface  and  redistributes 
the  momenta  of  the  slave  and  master  nodes  (411  cards) 


SLIDE2  This  subroutine  computes  the  cross-sectional  area  and  volume 
change  to  the  slave  elements  caused  by  the  master  surface  in- 
truding into  the  elements  (153  cards) 

START  This  subroutine  reads  the  initial  dynamic  or  detonation  condi- 

tions and  sets  the  variables  equal  to  the  conditions  which  exist 
at  impact  and/or  explosive  detonation  (145  cards) 

STRESS  This  subroutine  computes  stresses  and  energies  for  the  solid- 
liquid  elements.  The  stresses  consist  of  elastic,  plastic  and 
viscous  stresses,  hydrostatic  pressure  and  artificial  viscosity 
(155  cards) 

The  plotting  program  for  state  plots  reads  data  from  Tape  2 and  generates 
plots  for  geometry,  stress,  strain,  pressure  and  velocity  fields  at  specified 
times.  The  plotting  program  for  time  plots  reads  data  from  Tape  3 and  gen- 
erates time -dependent  plots  for  the  following  system  parameters:  center  of 
gravity,  kinetic  energy,  internal  energy,  total  energy,  plastic  work,  axial 
momentum,  axial  velocity,  spin  momentum,  spin  velocity,  maximum  axial 
coordinate  and  minimum  axial  coordinate.  Each  plot  presents  data  for  the 
projectile,  target  and  total  system. 


3.  2 NODE  AND  ELEMENT  ARRAYS 

There  are  13  node  variables  which  are  contained  in  common  arrays  and  19 
element  variables  contained  in  common  arravs.  The  node  arrays  are  de- 
scribed as  follows: 
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z 


T 

RDOT 

ZDOT 

TDOT 

FR 

FZ 

FT 

NODE 

AMASS 

PMASS 

IRZT 


The  radial  R coordinate  of  the  node  (see  Figure  2) 

The  axial  Z coordinate  of  the  node 

The  angular  6 coordinate  of  the  node 

The  radial  velocity  of  the  node 

The  axial  velocity  of  the  node 

The  angular  velocity  of  the  node 

The  radial  force  acting  on  the  node 

The  axial  force  acting  on  the  node 

The  circumferential  force  acting  on  the  node 

Identifies  the  node  number  such  that  it  is  not  necessary  to  have 
a consecutive  node  numbering  system.  The  node  number  must 
not  be  greater  than  the  dimension  of  the  node  arrays,  however. 

The  total  mass  of  the  node 

The  mass  of  the  node  contributed  by  the  projectile 

Identifies  radial,  axial  and  circumferential  restraints  on  the 
node.  It  also  identifies  a slave  node  and  indicates  if  it  is  free 
or  sliding  on  the  master  surface 


The  element  arrays  are  described  as  follows: 

NEL  Identifies  the  element  number.  Must  not  be  greater  than  the 

dimension  of  the  element  arrays. 


in 
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NODE1 

NODE  2 

NODE  3 

ICHECK 

SR 

SZ 

ST 

SRZ 

SRT 

SZT 

VOL 

DVOL 

MAT 

EBAR 

ES 


darea 


DDVOL 

TSTART 





Node  number  of  node  i (see  Figure  2) 

Node  number  of  node  j 
Node  number  of  node  m 

Identifies  state  of  material  (elastic,  plastic,  fractured) 

Total  normal  stress  in  the  radial  direction 
Total  normal  stress  in  the  axial  direction 
Total  normal  stress  in  the  circumferential  direction 
Shear  stress  in  the  R-Z  plane 
Shear  stress  in  the  R-8  plane 
Shear  stress  in  the  Z-0  plane 
Initial  volume  of  the  element 
Volumetric  strain 
Material  number 

Equivalent  plastic  strain  defined  in  Equation  (13) 

' 

Total  internal  energy  per  original  unit  volume  defined  in  Equa- 
tion (27) 

The  cross-sectional  area  of  the  master  surface  intruding  into 
a slave  element 

The  volume  of  the  master  surface  intruding  into  a slave  element 

The  time  required  for  the  detonation  wave  to  arrive  at  the  near- 
est node  of  the  element  from  the  point  of  detonation. 
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SECTION  4 

PROGRAM  USER  INSTRUCTIONS 


The  required  EPIC -2  input  data  for  various  conditions  are  summarized  in 
Figure  11  and  described  in  subsections  4.  1 through  4.  4.  Output  data  are 
described  in  subsection  4.  5.  Succeeding  subsections  deal  with  diagnostics, 
estimated  central  processor  time  requirements,  and  central  memory  stor- 
age requirements  and  alterations. 


4.  1 INPUT  DATA  FOR  AN  INITIAL  RUN 

An  initial  run  is  one  which  requires  input  data  for  the  initial  geometry  and 
impact/ detonation  conditions  at  time  = 0.  The  descriptions  which  follow 
are  for  the  input  data  in  Figure  11.  Consistent  units  must  be  used. 


Identification  Card  (415,  F10.0) 

CASE  = Case  number  for  run  identification 
CYCLE  ■ 0 for  initial  run 

IMAT  = A code  describing  the  types  of  materials.  Two  options  are 
available. 

IMAT  = 1 reads  material  properties  for  solid-liquid 
elements,  (no  explosive  element  data  read) 

IMAT  = 2 reads  material  properties  for  solid-liquid 
elements  and  explosive  elements. 
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INITIAL  RUN  INPUT  DATA 

IDENTIFICATION  CARD  1415.  F 10. 01 
CASE  |CYCU|  I MAT  jGEOM|  CPMAX 


23  SOL  ID -LIQUID  MATERIAL  CARDS  IX IS. 81 


! MATH 

MATL2 

MAIL) 

MATL4 

MATL5  { < 

( EXPLOSIVE  MATERIAL  CARDS  FOR  IMAT  • 2 (54 15  81 
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MATL10  ! 1 
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frlSLID[  IRIC[ 
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BLANK  CARD  $ 
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TARGET  NODE  DATA 
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BLANK  CARD 


PROJECTILE  ELEMENT  DATA  CARDS  ' AS  REQUIRED  1915) 

TOFT 


ELE1  | ELEN  | ^ 


MATL 

M 

M 

M 

BLANK  CARD 


PROJECTILE  ELEMENT  SPECIAL  SHAPES  - AS  REQUIRED  (DESCRIPTION  FOLLOWSI 
| BLANK  CARD  ( 

TARGET  ElfMENT  DATA  CARDS  - AS  REQUIRED  19151 


MATE  | N1 

N2  | N3 

M 

[TOq 

UmcJ 

BLANK 


CARD  } 


Figure  11.  Summary  of  Input  Data 


, 1 
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TARGET  ElEMENT  SPECIAL  SHAPES  - AS  REQUIRED  (DESCRIPTION  FOLLOWS) 
| BLANK  CARD  } 

CONCENTRATED  MASS  CARDS  - AS  REQUIRED  <15.  E15.8I 


| N | MASS  IN) 

RIGID  BODY 

IDENTIFICATION  CARDS 

- AS 
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(1515) 
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QdT 

iH 

SLIDING  SURFACE 
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REQUIRED  0151 
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NSE 
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LJ 
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rn 

1 S 15 
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1 S£16| 

INITIAL  VELOCITY /DETONATION 

CARD 

(8F10.0) 

| PRDOT 

| PZDOT 

| PTDOT 

TRDOT 

| TZDOT 

j nooT 

| RDET  1 

ZDET 

INTEGRATION 

TIME 

INCREMENT 

CARD 

(6F10.0) 

DTI 

| DTMAX 

| DTMIN 

SSF 

TMAX 

TBURN 

DATA  OUTPUT  CARDS  - 

AS  REQUIRED  (4F10.0,  3151 

| TIME 

ECHECK 

BCHECK 

RDAMP 

ILOAD 

ISAVE 

idataI 

RESTART  RUN  INPUT  DATA 


IDENTIFICATION  CARD  (315.  5X.  F10.0) 


[CASE  CYClEjlWINO^^ 

CPMAX 

INTEGRATION  TIME  INCREMENT  CARD  I3F10.0I 

| DTMIN  | SSF 

DATA  OUTPUT  CARDS  - AS  REQUIRED  (4F10.0,  315) 


TIME 


ECHECK 
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RDAMP  ILOAD  ll  SAVE!  IDATi 


Figure  11.  Summary  of  Input  Data  (Contirued) 
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SPECIAL  SHARES  INPUT  DATA  FOR  NODES 


||Q^3QQ| 


IKNTIF ICATION  CARD  FOR  ROD  NODES  141*.  *10.01 


ZTOP  ZBOT  EXPAND 


FOP  RADII  CARD  IS  I FOR  ROD  NODES  I8F10. 0) 


RTIll  RTI2I 


BOTTOM  RADII  CARO  (SI  FOR  ROD  NODES  (OF  10. 01 


RBIII  I RBl?l 


IDENTIFICATION  CARD  FOR  NOSE  NODES  HI*.  *X.  110.01 


/TOP 


TOP  RADII  CARDlSI  FOR  NOSE  NODES  I8F10.0I 


RTIll  I RT(?> 


MINIMUM  Z COORDINATE  CARDlSI  FOR  NOSE  NODES  I8F10.0I 


/Mill  I /mi/  I 


CARD  FOR  SPHERE  NODES  HI*.  *X.  2F10.0I 


/TOP  ZBOT 


IHBEEI 


IDENTIFICATION  CARD  FOR  FEAT  PUTE  NODES  121*1 


DESCRIPTION  CARD  FOR  FLAT  PLATE  NODES  141*.  6F10.0I 


RMAX 

| RMIN 

/MAX 

! /MIN 

R -EXPAND 

2 -EXPAND 

SPECIAL  SHAPES  INPUT  DATA  FOR  ELEMENTS 

IDENTIFICATION  CARD  FOR  ROD  ELEMENTS  (61*1 


MATERIAL  CARO  FOR  ROD  ELEMENTS  11*151 


[dUEflj 


iHgjl 


IDENTIFICATION  CARD  FOR  NOSE  ELEMENTS  HI*.  1QX.  1*1 


MATERIAL  CARD  FOR  NOSE  ELEMENTS  (1*1*1 


irmrqii 


CARD  FOR  SPHERE  ELEMENTS  HI*.  UK.  21*1 


CARD  FOR  FIAT  PLATE  ELEMENTS  171*1 


Figure  11.  Summary  of  Input  Data  (Continued) 
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STATE  PLOTS  INPUT  DATA 

STATE  PLOT  IDE NTIFICATION  CAROS  AS  REQUIRED  I? IS.  4F10.0) 


TYPE 


|CYCLt[ 


2MAX 


2MIN 


RMAX 


RMIN 


DEFORMED  GEOMETRY  PLOT  CARO  FOR  TYPE  • I (3A6.  EX,  4151 


TITLE 

■« 

Dl 

ua 

o 

STRESS  FIELD  PLOT  CARD  FOR  TYPE  • 2 

I3A6,  EX.  215.  SF10.0I 

TITU 

■« 

STRESSIll 

STRESSIEI 

STRESSUI 

STRESSI4I 

STRESSI5I 

PRESSURE  FIELD  1 

PLOT  CARD 

FOR  TYPE  • 

) I3A6.  ?X,  EI5,  5F10.0I 

TITLE 

■IBB 

imi 

PRESI1I 

PRESIEI 

PRESI3I 

PRESI4) 

PRESI5I 

STRAIN  FIELD  PLOT  CARD  FOR  TYPE  • 4 

I3A6,  EX.  EI5.  5F10.0I 

TITU 

i 

■ 

mi 

STRAIN(l) 

strainiei 

STRAINI3I 

STRAINI4I 

STRAINI5I 

VELOCITY  FIELD  PLOT  CARD 

FOR  TYPE  • 

5 I3A6,  EX,  15.  5X.  F10.0I 
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TOTAL  E FARCY 
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IR  VELOCITY  FOR  PLANE  STRAINI 


BLANK  CARD 


ENDS  RUN 


Figure  11.  Summary  of  Input  Data  (Concluded) 
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IGEOM 


= A code  describing  the  geometry  of  the  problem 
IGEOM  = 1 defines  axi-symmetric  geometry 


C PM  AX  = 


IGEOM  = 2 defines  plane  strain  geometry 

Central  processor  time  (hours)  at  which  the  results  will  be 
written  on  the  restart  and  plot  tape  (if  previously  requested) 
and  the  run  will  stop.  This  feature  will  be  bypassed  if 
CPMAX  = 0. 


23  Material  Cards  for  Solids  and  Liquids  (5E15.  8) 

These  cards  are  included  for  IMAT  - 1 or  IMAT  = 2.  Each  of  the  23  mat- 
erial cards  provides  a specific  material  characteristic  for  five  materials. 
For  instance,  the  first  card  describes  the  density  of  the  materials.  The 
density  of  material  1 is  entered  in  columns  1-15,  the  density  of  material 
2 in  columns  16-30,  etc.  The  second  card  contains  the  specific  heats  for 
the  five  materials,  etc.  It  is  only  necessary  to  describe  the  materials 
used  for  the  analysis.  If  materials  1 and  3 are  used  for  a specific  analysis, 
all  the  material  1 data  must  be  specified  in  columns  1-15  and  all  the  ma- 
terial 3 data  in  columns  31-45.  Columns  16-30,  46-60,  and  61-75,  repre- 
senting materials  2,  4 and  5,  can  be  left  blank. 

Card  1 Density  (mass /volume) 

Card  2 Specific  heat  (work/mass/degree) 

Card  3 Shear  modulus  of  elasticity  (force/area) 

Card  4 Absolute  viscosity  for  Navier-Stokes  Equations  (force-time/area) 
Card  5 Yield  stress  (force/area) 

Card  6 Ultimate  stress  (force/ area) 
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Card  7 


Card  8 
Card  9 
Card  10 
Card  11 
Card  12 
Card  13 
Card  14 
Card  15 
Card  16 
Card  17 
Card  18 
Card  19 
Card  20 
Card  21 

Card  22 

Card  23 


Strain  at  which  the  ultimate  stress  is  achieved.  Should  be  consist- 
ent with  the  equivalent  plastic  strain  definition  in  Equation  (13) 
which  is  essentially  a "true"  strain.  Stress  varies  linearly  between 
the  yield  stress  and  the  ultimate  stress.  The  stress  for  larger 
strains  is  equal  to  the  ultimate  stress  until  fracture  occurs. 

Maximum  negative  hydrostatic  pressure  (force/ area) 

Strain  rate  effect  constant,  Cj,  in  Equation  (24) 

Pressure  effect  constant,  C2>  in  Equation  (24) 

Pressure  effect  constant,  C 3>  in  Equation  (24) 

Temperature  effect  constant,  C4>  in  Equation  (24) 

Temperature  effect  constant,  C5,  in  Equation  (24) 

Hydrostatic  pressure  constant,  Kj,  in  Equation  (26)  (force/area) 

Hydrostatic  pressure  constant,  K2«  in  Equation  (26)  (force/area) 

Hydrostatic  pressure  constant,  K3>  in  Equation  (26)  (force/area) 

Gruneisen  coefficient,  r,  in  Equation  (26) 

Linear  artificial  viscosity  coefficient,  , in  Equation  (31) 

Quadratic  artificial  viscosity  coefficient,  , in  Equation  (31) 

Equivalent  strain  | If  either  is  exceeded  there  is  a shear  and 

Volumetric  strain  | tensile  failure-  Positive  hydrostatic  pres- 
sure and  viscosity  stress  capability  remain. 

If  Equivalent  strain  is  negative;  material 
behaves  like  a liquid. 

Equivalent  strain  If  exceeded  the  element  fails  totally  and 

produces  no  stresses  or  pressures. 

Initial  temperature  (degrees) 
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Six  Material  Cards  for  Explosives  (5E15.  8) 
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These  cards  are  included  only  if  IMAT  = 2 in  the  Identification  Card.  The 
explosive  material  numbers  are  designated  6 through  10  with  material  6 
entered  in  columns  1-15,  etc. 

Card  1 Density  (mass/volume) 

Card  2 Internal  Energy.  Initial  value  of  E in  Equation  (29)  (energy/volume) 
Card  3 Detonation  velocity  (distance/time) 

Card  4 Material  constant,  y , in  Equation  (29) 

Card  5 Linear  artificial  viscosity  coefficient.  C.  , in  Equation  (31) 

o 

Card  6 Quadratic  artificial  viscosity  coefficient.  Cq  , in  Equation  (31) 


Miscellaneous  Card  (715) 


NSLID 


TRIG 


NMASS 


NDEP 


= The  number  of  sliding  surfaces  (currently  limited  to  a max- 
imum of  five).  This  does  not  include  the  rigid  sliding  sur- 
face which  follows. 

= 1 gives  a rigid  frictionless  surface  on  the  positive  side  of  the 

plane  described  by  Z = 0.  If  the  equations  of  motion  cause  a 
node  to  have  a negative  Z coordinate,  the  node  coordinate  is 
redefined  and  set  to  Z = 0.  If  IR1G  = 0 this  option  is  not  used. 

= The  number  of  concentrated  masses  to  be  entered  late'r.  These 
masses  are  in  addition  to  the  concentrated  masses  automatic- 
ally generated  from  the  elements. 

= The  number  of  systems  of  nodes  which  are  designated  to  travel 
as  a rigid  body  in  the  axial  direction.  It  does  not  apply  to  the 
radial  velocities.  The  specific  nodes  are  input  later. 


l 
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The  next  three  variables  (NTOP,  NBOT,  NRING)  are  used  to  identify  the 
nodal  geometry  for  internal  load  calculations  in  a slender  projectile  with 
axi-symmetric  geometry  (IGEOM  = 1).  If  this  option  is  used,  the  nodal 
geometry  must  be  consistent  with  the  arrangement  of  nodes  generated  by 
the  special  shapes  generator  for  rod  geometry.  Leave  blank  if  it  will  not 
be  used, 

NTOP  = The  node  number  of  the  centerline  node  at  the  top  free  end 
of  the  cylindrical  projectile. 

NBOT  = The  node  number  of  the  centerline  node  at  the  lower  end  of 
the  cylindrical  portion  of  the  projectile. 

NRING  = The  number  of  rings  of  nodes  in  the  cylindrical  projectile. 

This  will  be  described  in  the  discussion  of  the  special  shapes. 


Projectile  Scale/Shift/Rotate  Card  (5F10.  0) 


RSCALE  = 


ZSCALE  = 
ROTATE  = 


RSHIFT 


Factor  by  which  the  R coordinates  of  all  projectile  nodes  are 
multiplied.  Applied  after  the  coordinate  shifts  (RSHIFT, 
ZSHIFT)  and  before  the  rotation  (ROTATE)  described  later. 

Factor  by  whiich  the  Z coordinates  are  multiplied. 

Rotation  about  R = Z = 0 in  the  R-Z  plane  of  all  projectile 
nodes  (degrees).  (Positive  clockwise)  Applied  after  the 
coordinate  shifts  (RSHIFT,  ZSHIFT),  and  the  scale  factors 
(XSCALE,  ZSCALE). 

Increment  added  to  the  R coordinates  of  all  projectile  nodes 
(length).  Applied  before  the  scale  factors  (RSCALE,  ZSCALE). 


ZSHIFT  = Increment  added  to  the  Z coordinates  (length). 


Projectile  Node  Data  Cards  (315,  2X,  311,  5F10.  0) 


The  Projectile  Node  Data  Cards  are  provided  as  required  for  the  projectile 
nodes.  If  a node  is  at  the  interface  of  the  projectile  and  the  target  and  con- 
tains mass  from  both  the  projectile  and  the  target,  it  must  be  included  with 
the  projectile  nodes.  The  node  numbers  (N1  . . . NN)  must  not  exceed 
the  dimension  of  the  node  arrays,  and  they  need  not  be  numbered  consecu- 
tively or  in  increasing  order.  End  with  a blank  card. 

N1  = First  node  in  a row  of  nodes. 

NN  = East  node  in  a row  of  nodes.  Leave  blank  if  a single  node 

(not  a row)  is  entered. 

INC  = Node  number  increment  between  adjacent  nodes.  If  left 

blank,  program  will  set  INC  = 1. 

IR  =1  will  restrain  all  nodes  (N2  • • • NN)  in  the  radial  R direc- 

tion. No  restraint  if  left  blank. 

1Z  =1  will  restrain  nodes  in  the  axial  Z direction 

IT  =1  will  restrain  nodes  in  the  circumferential  9 direction 

R1  = Radial  coordinate  of  node  N1  (length) 

% 

Z1  = Axial  coordinate  of  node  N1  (length) 

RN  = Radia  coordinate  of  node  NN  (length).  Leave  blank  if  a single 

node  is  entered. 

ZN  - Axial  coordinate  of  node  NN  (length) 

EXPAND  = Factor  by  which  the  distance  between  adjacent  nodes  changes. 

For  example,  if  the  distance  between  the  first  two  nodes  is 
Aj,  the  distance  between  the  second  and  third  nodes  is  A0  = 
Aj.  EXPAND.  A summary  of  relative  spacing  between  the 
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first  two  nodes  and  the  last  two  nodes  is  given  in  Figure  12. 
Analytic  expressions  are  also  provided  for  cases  not  included 
in  the  Figure.  If  left  blank,  the  program  will  set  EXPAND  = 
1.0  for  uniform  spacing. 


NODES  Nl  (Ml*  INC) 


(MM  - INC) 


INCREMENTS  A | 


TOTAL  LENGTH  _ 
USER  INCREMENTS  M 


A,.,  - A.  . EXPAND 


A1  N(l-EXPAND) 
A (1-EXPANDN) 


•VM  _N(1- 


J V N, 


11 ' (Txf>ANb>  > 

Figure  12.  Nodal  Spacing  for  Various  Expansion  Factors 

Projectile  Node  Special  Shapes 


These  are  described  later.  End  with  a blank  card. 
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Target/Scale/Shift/Rotate  Card  (5F10.  0) 

Same  as  Projectile  Scale/ Shift/ Rotate  card  except  it  applies  to  the  target 
nodes. 
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Target  Node  Data  Cards  (315,  2X,  311,  5F10.  0) 

Same  as  the  Projectile  Note  Data  cards.  If  there  are  interface  nodes  which 
contain  mass  from  both  the  projectile  and  the  target,  these  nodes  are  in- 
cluded with  the  projectile  nodes  and  must  not  be  included  here.  End  with  a 
blank  card. 


Target  Node  Special  Shapes 


These  are  defined  later.  End  with  a blank  card. 


Projectile  Element  Data  Cards  (915) 

The  Projectile  Element  Data  cards  are  provided  as  required  for  the  pro- 
jectile elements.  The  element  numbers  (ELE1  • • • ELEN)  must  not  ex- 
ceed the  dimension  of  the  element  arrays,  and  they  need  not  be  numbered 
consecutively.  End  with  a blank  card. 

ELE1  = First  individual  element  in  a series  of  elements. 

ELEN  = Last  individual  element  in  a series  of  elements.  Must  be 

greater  than  ELE1.  Leave  blank  if  a single  element  is 
entered. 
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= Element  number  increment  between  successive  elements. 

If  left  blank,  the  program  will  set  ELE  INC  = 1.  (Example: 
If  EL  El  = 100,  ELEN  = 104,  and  ELE  INC  = 2,  there  will 
be  three  elements  designated  at  100,  102,  and  104). 

MATL  = Material  number  (1-10)  of  the  elements.  If  left  blank,  the 
material  number  from  the  previous  element  data  card  will 
be  used. 


N1  - N4  = Node  numbers  for  the  first  composite  element  which  contains 
two  individual  triangular  elements.  Nodes  Nl,  N2  and  N3 
(Counterclockwise)  define  the  first  element  in  the  composite 
elements  and  nodes  Nl,  N3  and  N4  (counterclockwise)  define 
the  second  element.  If  N4  is  blank,  the  second  element  will 
not  be  generated. 

NODE  INC  = The  node  number  increment  added  to  the  node  numbers  of 

the  previous  composite  or  individual  element.  If  left  blank, 
the  program  will  set  NODE  INC  = 1.  (Example:  If  Nl  - N4 
are  24,  20,  22  and  26  and  if  NODE  INC  = 4,  then  the  nodes 
of  the  second  composite  element  are  28,  24,  26,  and  30). 


Projectile  Element  Special  Shapes 


These  are  defined  later.  End  with  a blank  card. 


Target  Element  Data  Cards  (915) 

Same  as  Projectile  Element  Data  cards.  End  with  a blank  card. 
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Target  Element  Special  Shapes 

These  are  defined  later.  End  with  a blank  card. 


Concentrated  Mass  Cards  (15.  Elf).  a) 

There  are  NMASS  (Defined  in  Miscellaneous  card)  cards  entered  for  (he 
concentrated  masses.  Each  card  contains  data  for  one  mass. 

N = The  nodi'  number  to  which  the  concentrated  mass  is  added. 

MASS(N)  - The  concentrated  mass  added  to  nodi'  N. 

Rigid-Hodv  Identification  Cards  (15) 

Each  system  of  rigid-body  nodes  contains  one  Rigid-Body  Identification 
card  and  one  Rigid-Body  Nodes  card.  The  program  is  dimensioned  for  a 
maximum  of  five  rigid-body  systems  wit.li  a maximum  15  nodes  per  system. 

If  there  are  no  rigid-body  systems  (NI)EP  0 in  Miscellaneous  card)  this 
group  of  cards  is  omitted.  If  there  is  more  than  one  system  of  rigid-body 
nodes,  both  cards  for  the  system  are  entered  before  the  Rigid-Body  Iden- 
tification card  for  the  second  system  is  entered. 

] 

i 1 

NI)N  - The  number  of  rigid-body  nodes  in  this  system. 

Rigid- Body  Nodes  Card  (1515)  ] 

11)1,  ID2,  - . . are  the  nodi'  numbers  for  the  rigid-body  nodes  which  can  be 
input  in  any  order.  These  nodes  will  always  have  identical  velocities  in  the 
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Z direction.  Master  nodes  (defined  later)  can  be  included  but  slave  nodes 
(defined  later)  cannot  be  included. 


Sliding  Surface  Identification  Cards  (315) 

Each  sliding  surface  contains  one  identification  card  and  cards  (as  required) 
describing  the  master  nodes,  slave  nodes  and  the  slave  elements.  The  pro- 
gram is  dimensioned  for  a maximum  of  five  sliding  surfaces  which  can  each 
contain  a maximum  of  100  master  nodes,  100  slave  nodes  and  100  slave 
elements.  If  there  are  no  sliding  surfaces  (NSLID  = 0 on  Miscellaneous 
card)  this  group  of  cards  is  omitted.  If  there  is  more  than  one  sliding 
surface,  all  data  for  the  first  sliding  surface  are  entered  before  the  second 
Sliding  Surface  Identification  card  is  entered. 

NMN  = The  number  of  master  nodes  on  the  sliding  surface 

NSN  = The  number  of  slave  nodes  on  the  sliding  surface 

NSE  = The  number  of  slave  elements  on  the  sliding  surface.  A slave 

element  is  a triangular  element  which  must  contain  two  slave  nodes. 

If  an  element  is  designated  as  a slave  element,  the  cross-sectional 
area  (and  the  resulting  volume)  will  be  adjusted  if  a master  node  is 
located  between  the  two  slave  nodes.  It  is  not  necessary,  however, 
to  include  slave  elements  to  allow  the  problem  to  run. 

I 

.Master  Node  Cards  (1515) 


IM1,  IM2,  . . . are  the  node  numbers  of  a row  of  master  nodes.  The  nodes 
must  be  entered  in  order  from  the  first  master  node  to  the  last  master  node 
along  the  row  of  nodes.  When  moving  from  the  first  node  to  the  last  node, 
the  slave  nodes  (and  elements)  must  be  to  the  left  of  the  master  surface. 
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Slave  Nodi'  Cards  (1615) 


IS1,  IS2,  . . . are  the  node  numbers  of  the  slave  nodes.  No  special  order 
is  required.  If  there  are  no  slave  elements  (NSE  0)  it  is  possible  to  desig 
nate  interior  nodes  (as  well  as  surface  nodes)  as  slave  nodes.  This  proce- 
dure allows  elements  containing  slave  nodes  to  fail  completely  to  simulate 
the  eroding  process.  The  slave  node  spacing  should  not  be  significantly 
greater  than  the  master  node  spacing  or  localized  deformations  in  the  mas- 
ter surface  may  occur.  A slave  node  can  be  restrained  only  in  the  radial 
direction,  and  only  if  the  radial  coordinate  of  the  slave  node  is  equal  to  the 
radial  coordinate  of  the  first  or  last  master  node,  which  must  also  be  re- 
strained in  the  radial  direction.  This  allows  a slavo  node  on  the  Y,  axis  to 
be  restrained. 


Slave  Element  Cards  (1615) 


ISE1  . 1SE2,  . . . are  the  element  numbers  of  the  slave  elements.  If  very 
severe  distortions  are  expected  on  the  slave  surface,  it  may  be  desirable 
to  not  designate  anv  slave  elements  since  numerical  instability  problems 
may  occur.  If  an  element  is  designated  as  a slave  element  , the  element 
must  contain  exactly  two  nodes  which  are  slave  nodes. 


Initial- Velocity/ detonation  Card  (tlKlO  0) 

This  card  describes  the  initial-velocity  and/or  detonation  conditions.  If 
there  are  interface  nodes  which  include  mass  from  both  the  projectile  and 
thi'  target,  the  velocities  of  these  nodes  are  adjusted  by  the  program  to 
conserve  momenta. 


PRDOT 

PZDOT 

PTDOT 

TRDOT 

TZDOT 

TTDOT 

RDET 

ZDET 


Projectile  velocity  in  the  R direction  (distance/time) 

Projectile  velocity  in  the  Z direction 

Projectile  velocity  in  the  0 direction  (radians/time) 

Target  velocity  in  the  R direction 

Target  velocity  in  the  Z direction 

Target  velocity  in  the  0 direction 

Radial  coordinate  of  the  detonation  (distance) 

Axial  coordinate  of  the  detonation 


Integration  Time  Increment  Card  (6F10.  0) 

DTI  = Integration  time  increment  (At  in  equations  of  motion  in  sub- 

section 2.  5)  for  the  first  cycle.  This  must  be  less  than  the 
time  required  to  travel  across  the  minimum  altitude  of  each 
triangular  element  at  the  sound  speed  of  the  material  in  that 
element. 

DTMAX  = The  maximum  integration  time  increment  which  will  be  used 

for  the  equations  of  motion.  If  At  from  Equation  (42)  is  greater 
than  DTMAX.  it  will  be  redefined  at  At  = DTMAX. 

DTMIN  = The  minimum  integration  time  increment  allowed.  If  At  from 
Equation  (42)  is  less  than  DTMIN,  the  results  will  be  written 
onto  the  restart  tape  and  the  run  will  stop. 

SSF  = The  fraction  of  the  sound  speed  transit  time  used  for  the  inte- 

gration time  increment.  This  is  identical  to  Ct  in  Equation 
(42)  and  must  be  less  than  unity. 
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TMAX  = The  maximum  time  the  problem  is  allowed  to  run.  This  time 
refers  to  the  dynamic  response  of  the  system,  not  the  central 
processor  time  (CPMAX)  defined  in  the  Identification  card. 

TBURN  = The  time  at  which  the  detonation  begins  at  RDET,  ZDET. 

Data  Output  Cards  (4F10.  0,  315) 

These  cards  are  used  to  specify  various  forms  of  output  data  at  selected 

times.  The  last  card  must  be  for  a time  greater  than  TMAX,  even  though 

output  will  not  be  provided  for  that  specific  time. 

TIME  = Time  at  which  output  data  will  be  provided. 

ECHECK  = A code  which  governs  the  printed  output.  The  three  following 
options  are  provided; 

1)  If  ECHECK  is  greater  than  1000.  , the  individual 
node  data  and  element  data  will  not  be  printed. 

Only  system  data  such  as  eg  positions,  momenta, 
kinetic  energies  and  average  velocities  are  pro- 
vided for  the  projectile,  target  and  entire  system. 

2)  If  ECHECK  is  positive  (includes  0.  ) and  less  than 
1000.  , the  system  data  and  individual  node  data 
will  be  printed.  Individual  solid-liquid  element 
data  will  be  printed  for  all  elements  which  have 
an  equivalent  plastic  strain  (Equation  (13))  equal 
to  or  greater  than  ECHECK.  For  example,  if 
ECHECK  = 0.  , all  element  data  will  be  printed. 

If  ECHECK  = 0.  5,  only  those  elements  with 
equivalent  plastic  strains  equal  to  or  greater 


BCHECK 


RDAMP 


I LOAD 


ISAVE 

I DATA 


than  0.  5 will  have  data  printed.  If  ECHECK  = 999.  , 
no  element  data  will  be  printed. 

3)  If  ECHECK  is  negative,  the  system  data  and  indiv- 
idual node  data  will  be  printed.  The  individual 
element  data  will  be  printed  only  for  those  elements 
which  are  plastic  or  failed.  Data  for  the  elastic 
elements  will  not  be  printed. 

A code  which  governs  the  printed  output  for  explosive  elements. 
Individual  element  data  will  be  printed  for  all  explosive  ele- 
ments which  have  a burn  fraction  (Equation  (30))  equal  to  or 
greater  than  BCHECK.  If  BCHECK  = 0.  all  explosive  elements 
data  will  be  printed.  If  BCHECK  = 0.  5,  only  those  explosive 
elements  with  burn  fractions  equal  to  or  greater  than  0.  5 will 
have  data  printed. 

The  radial  damping  constant,  C^,  in  Equation  (38).  This  damp- 
ing term  acts  until  the  time  specified  in  the  following  Data  Out- 
put card. 

1 will  compute  internal  loads  in  a cylindrical  projectile.  This 
option  can  only  be  used  if  NTOP,  NBOT  and  NRING  are  properly 
defined  in  the  Miscellaneous  card.  Will  not  compute  if  left 
blank. 

1 will  write  results  on  Tape  2 for  possible  restart  runs  or  state 
plots.  Will  not  write  results  if  left  blank. 

1 will  write  system  data  on  Tape  3 for  possible  time  plots. 

Will  not  write  system  data  if  left  blank. 


With  the  exception  of  the  special  shapes,  this  completes  the  description  of 
the  input  data  for  an  initial  run.  Various  types  of  special  shapes  can  be 
generated  for  either  the  projectile  or  the  target.  The  node  and  element 
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data  for  these  special  shapes  are  entered  individually  in  the  locations  iden 
tified  in  Figure  11.  A description  of  the  special  shapes  data  follows. 


i 


Identification  Card  for  Rod  Nodes  (415,  3F10.  0) 

The  rod  geometry  descriptions  are  given  in  Figure  13.  The  axi-symmetric 
rod  geometry  is  generated  about  the  Z axis  as  shown.  If  the  inner  radius 
coincides  with  the  Z axis  (R  = 0),  a solid  rod  is  generated  and  the  radius 
is  restrained  in  the  R direction.  If  the  inner  radius  does  not  coincide  with 
the  Z axis,  a hollow  rod  is  generated  and  there  are  no  restraints. 


1 

N1 


NRING 

NPLN 
ZTOP 
ZBOT 
EX PAN  D 


Identification  number  for  rod  geometry. 

The  node  number  of  the  inner  radius  node  at  the  top  surface 
of  the  rod.  The  nodes  of  the  rod  are  numbered  consecutively 
beginning  with  Nl.  They  are  numbered  in  horizontal  rows  (in- 
creasing in  the  radial  direction)  beginning  at  the  top  and  work- 
ing downward. 

The  number  of  vertical  rings  of  nodes.  The  maximum  allow- 
able number  of  rings  is  16. 

The  number  of  horizontal  rows  of  nodes. 

The  Z coordinate  at  the  top  of  the  rod. 

The  Z coordinate  at  the  bottom  of  the  rod. 

Factor  by  which  the  vertical  distance  between  adjacent  nodes 
changes.  Factor  applies  from  top  to  bottom.  Same  as  de- 
scribed for  the  Projectile  Node  Data  cards. 
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Bottom  Radii  Cards  for  Rod  Nodes  (8F10.  0) 

r 

RB  (1)  - RB  (NRING)  = Radii  of  the  rings  at  the  bottom  of  the  rod. 

Identification  Card  for  Nose  Nodes  (315,  5X,  F10.  0) 

The  nose  geometry  descriptions  are  given  in  Figure  14.  The  axi-symmetric 
geometries  for  rounded  and  conical  nose  shapes  are  also  generated  about 
the  Z axis  (pointed  downward)  and  the  centerline  nodes  are  restrained  in 
the  radial  direction.  The  nodes  at  the  rod-nose  interface  are  not  generated 
with  the  nose  generator  and  must  therefore  be  generated  with  the  rod 

generator. 
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ZTOP 
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!2  Identifies  a rounded  nose  shape. 

3 Identifies  a conical  nose  shape. 

= The  node  number  of  the  centerline  node  at  the  top  of  the  nose 
at  the  rod  interface.  This  is  identical  to  the  bottom  center- 
line  node  of  the  rod.  The  nodes  are  numbered  consecutively 
in  a clockwise  direction  (excluding  the  interface  nodes)  begin- 
ning with  the  entire  center  ring  and  expanding  outward. 

= Number  of  rings  of  nodes  in  the  nose  shape.  Note  that  the 
number  of  rings  of  nodes  in  the  nose  is  one  less  than  the  num- 
ber of  rings  of  nodes  in  the  matching  rod  at  the  interface. 

The  maximum  allowable  number  of  rings  is  15. 

= The  Z coordinate  at  the  top  of  the  nose.  This  is  identical  to 
ZBOT  for  the  rod  shape  at  the  rod  interface. 


ROUNDED  NOSE  (ISHAPE-2) 
y-ROD 


[N1  + NRING  + 1] 
ZTOP  — 
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RT(1>  RTI2) 
1 2 


RT(NRING) 

I 

NRING 
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NODES 


NODES  AND  ELEMENTS  ARE  NUMBERED  CONSECUTIVELY,  IN  A 
CLOCKWISE  DIRECTION,  BEGINNING  WITH  THE  ENTIRE  CENTER 
RING  AND  EXPANDING  OUTWARD.  THIS  DOES  NOT  APPLY  TO 
THE  INTERFACE  NODES  SINCE  THEY  ARE  NOT  GENERATED  FOR 
THE  NOSE  SHAPES,  BUT  ARE  USED  ONLY  FOR  REFERENCE. 

Figure  14.  Nose  Shape  Geometry 


Top  Radii  Cards  for  Nose  Nodes  (8F10.  0) 

RT  (1)  - RT  (NRING)  = Radii  of  the  rings  at  the  top  of  the  nose.  Only  solid 
(not  hollow)  nose  shapes  are  allowed.  Note  that  the  centerline  node  is  al- 
ways on  the  Z axis  and  that  RT  (!)  is  the  adjacent  node. 


Minimum  Z Coordinate  Cards  for  Nose  Nodes  (8F10.  0) 

ZM  (1)  - ZM  (NRING)  = Minimum  Z coordinates  (at  R = 0)  for  the  rings  in 
the  nose  shape. 


Card  for  Sphere  Nodes  (315,  5X,  2F10.  0) 

The  sphere  geometry  descriptions  are  given  in  Figure  15.  The  axi-symmetric 
sphere  geometry  is  generated  about  the  Z axis  and  the  centerline  nodes  are  re- 
strained in  the  radial  direction. 


4 

N 1 


NRING 

ZTOP 

ZROT 


= Identification  number  for  sphere  geometry. 

= The  node  number  at  the  center  of  the  sphere.  The  nodes  are 
numbered  consecutively,  in  a clockwise  direction,  beginning 
with  the  inner  ring  and  expanding  outward  with  equal  radial 
spacing. 

= Number  of  rings  of  nodes.  The  maximum  allowable  number  of 
rings  is  15. 

= The  Z coordinate  at  the  top  of  the  sphere 

= The  Z coordinate  at  the  bottom  of  the  sphere. 
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NODES  AND  ELEMENTS  ARE  NUMBERED 
CONSECUTIVELY,  IN  A CLOCKWISE 
DIRECTION,  BEGINNING  WITH  THE 
CENTER  RING  AND  EXPANDING  OUTWARD 


Figure  15.  Sphere  Geometry 


i 
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Identification  Card  for  Flat- Plate  Nodes  (215) 

The  flat-plate  geometry  descriptions  are  given  in  Figure  16.  The  flat-plate 
geometry  is  generated  for  a rectangular  cross  section  as  shown. 

5 = Identification  number  for  flat- plate  geometry. 
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N1  = The  node  number  of  the  first  node  at  R = RMIN  and  Z = ZMAX.  The 
nodes  are  generated  in  rows  parallel  to  the  R axis  and  are  numbered 
consecutively  within  each  row,  in  the  direction  of  the  increasing  R 
axis. 


Figure  16.  Flat-Plate  Geometry 


Description  Card  for  Flat-Plate  Nodes  (415,  6F10.  0) 

NR  = The  number  of  nodes  in  the  R direction. 

NZ  = The  number  of  nodes  in  the  Z direction. 
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IDZ 


IRFIX 

RMAX 

RMIN 

ZMAX 

ZMIN 

R-EXPAND 

Z-EXPAND 


The  node  number  increment  between  adjacent  nodes  when 
moving  downward  parallel  to  the  Z axis.  If  the  entire  plate 
is  formed  with  one  node  shape,  the  IDZ  should  be  equal  to 
NR.  For  multiple-plate  shapes  it  should  be  equal  to  the 
total  number  of  nodes  in  the  R direction. 

1 will  radially  restrain  all  nodes  on  the  Z axis  (R  = 0). 

The  maximum  coordinate  in  the  R direction. 

The  minimum  coordinate  in  the  R direction. 

The  maximum  coordinate  in  the  Z direction. 

The  minimum  coordinate  in  the  Z direction. 

Factor  by  which  the  distance  between  adjacent  nodes  changes 
in  the  increasing  R direction  from  RMIN  to  RMAX. 

Factor  by  which  the  distance  between  adjacent  nodes  changes 
in  the  decreasing  Z direction  from  ZMAX  to  ZMIN. 


Identification  Card  for  Rod  Elements  (615) 


1 

N1 

NRCOL 


NZLAY 

IDIA 


= Identification  number  for  rod  geometry. 

= The  node  number  of  the  inner  radius  node  at  the  top  surface 
of  the  rod. 


= The  number  of  vertical  columns  of  elements.  Note  that 
NRCOL  = NRING  - 1 for  rod  geometry.  The  maximum 
allowable  number  of  columns  of  elements  is  15. 

= The  number  of  horizontal  layers  of  elements.  Note  that 
NZLAY  = NPLN  - 1. 


= Code  for  orientation  of  diagonals  as  shown. 


1 
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ELE1 


= The  element  number  of  the  first  element  (contains  node  Nl). 
The  elements  are  numbered  consecutively  and  are  generated 
in  vertical  columns  beginning  with  the  top  layer  1 and  ending 
at  bottom  layer  NZLAY.  The  entire  inner  radial  column  is 
generated  before  the  adjacent  outer  radius  column  is  gener- 
ated etc. 


Material  Card  for  Rod  Elements  (1515) 


M (1)  - M (NRCOL)  = Material  numbers  for  the  columns  of  elements. 


Identification  Card  for  Nose  Elements  (315,  IPX,  15)  ' 

The  nose  geometry  descriptions  are  given  in  Figure  14. 

12  identifies  a rounded  nose  shape. 

3 identifies  a conical  nose  shape. 

Nl  = The  node  number  at  the  centerline  node  at  the  top  of  the 

nose  at  the  rod  interface.  This  is  identical  to  the  bottom 
centerline  node  of  the  rod. 

NRING  = Number  of  rings  of  elements  (and  nodes)  in  the  nose  shape. 

The  maximum  allowable  number  of  tings  is  15. 

ELE1  = The  element  number  of  the  first  element  in  the  nose  shape. 

The  elements  are  numbered  consecutively,  in  a clockwise 
direction,  beginning  with  the  entire  center  ring  and  expand- 
ing outward. 
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Material  Card  for  Nose  Elements  (1515) 

M (1)  - M (NR1NG)  = Material  numbers  for  the  rings  of  elements. 


Card  for  Sphere  Elements  (315,  IPX.  215) 

The  sphere  geometry  descriptions  are  given  in  Figure  15. 
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N1 

NR1NG 

EL  El 


MATL 


= Identification  number  for  sphere  geometry. 

= The  node  number  at  the  center  of  the  sphere. 

= Number  of  rings  of  elements.  The  maximum  allowable  num- 
ber of  rings  is  15. 

= The  element  number  of  the  first  element  in  the  sphere.  The 
elements  are  numbered  consecutively,  in  a clockwise  direc- 
tion, beginning  with  the  inner  ring  and  expanding  outward. 

= The  material  number  for  all  elements  in  the  sphere. 


Pard  for  Flat- Plate  Elements  (715) 


The  flat- plate  descriptions  are  given  in  Figure  16. 
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5 

N1 

NRLAY 

NZLAY 


= Identification  number  for  flat- plate  geometry. 

= The  node  number  of  the  first  node  at  R = RMIN  and  Z = ZMAX 
= Number  of  vertical  layers  (columns)  of  elements. 

= Number  of  horizontal  layers  of  elements. 
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MATE 


= Code  for  orientation  of  diagonals  as  shown. 

= The  element  number  of  the  first  element  in  the  plate  (contains 
node  Nl).  The  elements  are  numbered  consecutively  in  hori- 
zontal rows  expanding  outward.  After  a row  is  complete,  the 
next  lower  row  is  numbered  in  a similar  manner. 

= The  material  number  for  all  elements  in  the  plate. 


4.2  INPUT  DATA  FOR  A RESTART  RUN 

A restart  run  is  one  which  reads  data  from  Tape  2.  These  data  must  have 
been  previously  generated  by  setting  ISAVE  = 1 on  a Data  Output  card  for  an 
initial  geometry  run  or  a previous  restart  run.  The  descriptions  which 
follow  are  for  the  input  data  in  Figure  11. 


Identification  Card  (315,  5X,  F10.  0) 


CASE 


CYCLE 


(WIND  = 


= Case  number  for  run  identification.  This  does  not  have  to 
be  identical  to  that  of  the  previous  run. 

= The  cycle  number  at  which  the  restart  occurs.  The  cycle 
numbers  for  which  restart  files  are  written  are  given  in  the 
printed  output  of  the  previous  run. 

0 will  not  rewind  Tape  2 after  the  program  has  read  the 
restart  data.  This  allows  the  entire  run  (from  time  = 0) 
to  be  saved  on  one  tape. 

< 

1 will  rewind  Tape  2 after  the  program  has  read  the  restart 

data.  Previously  generated  data  for  lower  cycle  numbers 

will  then  be  erased, 

v. 
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CPMAX  = Central  processor  time  (hours)  at  which  the  results  will  be 
written  onto  a restart  tape  and  the  run  will  stop.  This  fea- 
ture will  be  by-passed  if  CPMAX  = 0. 


Integration  Time  Increment  Card  (3F10.  0) 


DTMIN 


TMAX 


The  minimum  integration  time  increment  which  is  allowed. 

If  At  from  Equation  (42)  is  less  than  DTMIN,  the  results  will 
be  written  onto  the  restart  tape  and  the  run  will  stop. 

The  fraction  of  the  sound  speed  transit  time  used  for  the  in- 
tegration time  increment.  This  is  identical  to  Ct  in  Equation 
(42)  and  must  be  less  than  unity. 

The  maximum  time  the  problem  is  allowed  to  run.  This  time 
refers  to  the  dynamic  response  of  the  system,  not  the  central 
processor  time  (CPMAX)  defined  in  the  Identification  Card. 


Data  Output  Cards  (4F10,  0,  315) 

Identical  to  the  Data  Output  cards  for  the  initial  run  input  data  described  in 
subsection  4.  1. 


4.  3 INPUT  DATA  FOR  STATE  PLOTS 

Input  data  for  state  plots  are  shown  in  Figure  11.  This  plotting  program 
reads  data  from  Tape  2 which  must  be  previously  generated  by  setting 
ISAVE  = 1 on  a Data  Output  Card  for  an  initial  geometry  run  or  a restart 
run.  Five  types  of  plots  are  available;  each  plot  requires  a State  Plot 
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Identification  card  and  another  card  to  input  the  data  for  the  specific  type 
of  plot  requested.  The  data  must  be  input  in  groups  of  two  cards  and  the 
requested  cycle  numbers  must  be  equal  to  or  greater  than  the  previously 
requested  cycle  number.  Program  ends  with  a blank  card. 


State  Plot  Identification  Card  (215,  4F10.  0) 


This  identifies  the  type  of  plot  and  the  coordinates. 


TYPE 


CYCLE 


/.MAX 


A code  to  specify  the  type  of  plot  requested. 

TYPE  = 1 gives  a geometry  plot 
TYPE  = 2 gives  a stress  field  plot 
TYPE  = 3 gives  a pressure  field  plot 
TYPE  = 4 gives  a strain  field  plot 
TYPE  = 5 gives  a velocity  field  plot 

The  cycle  number  of  the  plot  which  is  desired.  The  cycle  num- 
bers of  the  data  written  on  Tape  2 are  given  in  the  printed  output 
of  the  main  program. 

The  maximum  axial  Z coordinate  on  the  plot.  The  Z coordinate 
(ZMAX  to  ZMIN)  is  10.  0 inches  long  and  is  divided  into  10  equal 
increments. 


ZMIN  = The  minimum  Z coordinate  on  the  plot. 

RMAX  = The  maximum  radial  R coordinate  on  the  plot.  The  R coordinate 

(RMAX  to  RMIN)  can  be  any  length  and  the  scale  is  equal  to  that 
of  the  Z coordinate. 


= The  minimum  R coordinate  on  the  plot. 


RMIN 


Deformed  Geometry  Plot  Card  (3A6,  2X,  415) 


This  gives  plots  of  geometry  composed  of  triangular  elements. 


TITLE 

EDGE 


IE 

ID 

IF 


The  title  which  is  written  on  the  plot.  The  time  and  cycle  num- 
ber are  also  written  on  the  plot. 

A code  to  specify  if  an  outline  of  the  geometry  should  be  plotted. 
EDGE  = 0 gives  no  outline. 

EDGE  = 1 gives  an  outline  of  the  actual  geometry. 

EDGE  = -1  gives  an  outline  of  the  mirror  image  of  the  actual 
geometry  for  axi-symmetric  cases. 

EDGE  = 2 gives  an  outline  of  both  the  actual  geometry  and 
the  mirror  image  of  the  actual  geometry. 

1 will  write  "E"  at  the  center  of  all  elements  which  are  in  the 
elastic  range.  Will  not  write  if  it  is  left  blank. 

1 will  write  "P"  at  the  center  of  all  elements  which  are  in  the 
plastic  flow  range.  Will  not  write  if  it  is  left  blank. 

1 will  write  "F"  at  the  center  of  all  elements  which  are  failed 
(no  shear  or  tensile  stress).  Will  not  write  if  it  is  left  blank. 


Stress  Field  Plot  Card  (3A6,  2X,  215,  5F10.  0) 

This  gives  plots  of  lines  of  constant  equivalent  stress  given  by  Equation  (23). 
TITLE  = The  title  which  is  written  on  the  plot. 

EDGE  = A code  to  specify  if  the  outline  of  the  geometry  should  be  plotted. 
See  description  given  for  the  Deformed  Geometry  Plot  card. 
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4.4  INPUT  DATA  FOR  TIME  PLOTS 


Input  data  for  time  plots  are  shown  in  Figure  11.  This  plotting  program 
reads  data  from  Tape  3 which  must  be  previously  generated  by  setting 
IDATA  = 1 on  the  Data  Output  cards  for  an  initial  run  or  a restart  run. 
The  variables  plotted  are  shown  in  Figure  11.  They  are  plotted  as  a 
function  of  time  at  the  various  times  specified  on  the  Data  Output  cards. 
Any  or  all  of  the  11  varables  can  be  plotted.  Each  plot  contains  data  for 
the  projectile,  the  target,  and  the  total  system  (projectile  plus  target). 
End  with  a blank  card. 


TYPE 

TMAX 

TMIN 

VMAX 

VMIN 

TITLE 


A code  to  specify  the  variable  to  be  plotted.  Code  is  given  in 
Figure  11. 

The  maximum  time  included  in  the  plot  (micro-seconds). 

The  minimum  time  included  in  the  plot  (mico-seconds). 

The  maximum  value  of  the  specified  variable  to  be  plotted. 
The  minimum  value  of  the  specified  variable  to  be  plotted. 
The  title  which  is  printed  on  the  plot. 


4.  5 OUTPUT  DATA  DESCRIPTION 

The  output  data  are  generally  described  using  the  terminology  of 
this  report.  A summary  of  the  printed  output  for  an  initial  run  is  as 
follows: 

• All  input  data  are  printed  exactly  as  read. 

• Material  data  are  printed  for  the  input  data. 
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• All  the  geometric  input  data  are  printed  in  a form  similar 
to  that  used  to  read  the  data. 

• For  each  node,  the  initial  coordinates,  mass  and  restraints 
are  printed.  The  slave  nodes  are  identified  by  adding  a "1" 
before  the  RZT  restraint.  (Example:  If  the  RZT  restraint 
is  100,  the  nodes  are  restrained  in  the  radial  direction.  If 
the  RZT  restraint  is  1100,  it  is  a slave  node  with  the  same 
restraint.  ) 

• For  each  element,  the  three  nodes,  cross-sectional  area, 
volume  and  material  number  are  printed. 

• A summary  of  the  initial  geometry  and  impact/detonation 
conditions  is  printed.  Included  are  the  number  of  nodes 

s and  elements,  mass,  eg  positions,  velocities,  detonation 

i 

location,  energies  and  integration  time  increment  data. 

i 

• For  each  integration  time  increment  selected  data  are 
printed.  These  consist  of  cycle  number  (CYCLE),  the 
time  (TIME),  the  integration  time  increment  (ITT).  the 
element  number  which  governs  the  time  increment 
(ECRIT),  the  total  system  kinetic  energy  (TOTAL  K.  E.  ). 
the  projectile  kinetic  energy  (PRO.T  K.  E.  ).  the  total  axial 
Z momentum  of  the  system  (TOTAL  MYZ),  the  projec- 
tile axial  Z momentum  (PRO.T  MVZ).  the  total  spin  mo- 
mentum of  the  system  (TOTAL  MVS),  the  projectile  spin 
momentum  (PRO.T  MVS),  and  the  total  energy  in  the 
system  (TOT  ENERGY).  If  ECRIT  = 0.  the  integration 
time  increment  calculated  by  Equation  (42)  is  either 
greater  than  DTMAX,  or  greater  than  1.  1 times  the  inte- 
gration time  increment  for  the  previous  cycle.  If  the 
problem  is  plane  strain  (instead  of  axi-symmetric), 
momenta  in  the  R direction  are  given  in  place  of  spin 
momenta  (TOTAL  MVS,  PRO.T  MVS). 
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• Node  data  are  printed  at  the  selected  times  specified  on  the 
Data  Output  Cards.  These  data  include  the  node  number 
(NODE),  the  coordinates  (R.Z.T),  the  velocities  (RDOT. 

ZDOT,  TDOT),  and  the  accelerations  (RDD,  ZDD,  TDD). 

If  the  acceleration  data  are  not  printed,  this  indicates  the 
node  is  a slave  node  which  is  currently  in  contact  with  the 
master  surface. 

• Element  data  are  also  printed  at  the  selected  times  on  the 
Data  Output  cards.  These  data  consist  of  the  element  num- 
ber (ELE),  the  relative  volume  (RELVOL).  the  volumetric 
strain  (DVOI.).  the  volumetric  strain  rate  (DVDOT).  the 
equivalent  plastic  strain  (EPBAR)  given  in  Equation  (13). 
the  equivalent  stress  (SEFF)  given  in  Equation  (23).  the 
normal  deviator  stresses  (SR.  SZ,  ST)  given  in  Equations 
(17)  - (19).  the  shear  stresses  (SRZ,  SRT.  SZT)  given  by 
Equations  (20)  - (22).  the  pressure  (PRESURE)  given  by 
Equation  (26)  or  (29),  the  artificial  viscosity  [Q  (VIS)  ] 
given  by  Equation  (31).  the  temperature  (TEMP)  given  by 
Equation  (25)  and  the  burn  fraction  (BURNFR)  given  by 
Equation  (30). 

- 

• Internal  loads  can  be  obtained  for  axi-symmetric  rod  pro- 
jectiles at  selected  times  if  the  nodal  geometry  is  consistent 
with  the  geometry  generator  for  rods  and  if  II.OAD  = 1 on 
the  Data  Output  cards.  The  output  is  given  along  the  center- 
line  of  the  projectile  from  nodes  NTOP  to  NBOT  as  specified 
on  t'n  .,ell.m-  card  Included  are  the  node  numbers 
(N«'l>E),  the  axial  coordinate  of  the  node  (Z),  the  compressive 
axial  load  (AXIAL  LOAD)  and  the  torque  (TORQUE).  The  re- 
sults are  valid  only  along  the  free  portion  of  the  projectile  and 
are  not  valid  for  the  portions  which  are  in  contact  with  the 
target. 
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• A summary  of  system  data  is  printed  for  selected  time 
| specified  on  the  Data  Output  cards.  Data  are  given  for  the 

projectile,  target  and  total  system.  These  data  consist  of 
mass,  eg,  kinetic  energy,  internal  energy,  total  energy 
(kinetic  plus  internal),  plastic  work,  momenta  in  Z direction, 
net  velocity  in  Z direction,  angular  momenta  about  Z axis 
(linear  momenta  in  R direction  if  plane  strain),  angular 
velocity  about  Z axis  (linear  velocity  in  R direction  if  plane 
strain),  maximum  Z coordinate  and  minimum  Z coordinate. 

4.  (i  DIAGNOSTICS 

The  following  diagnostics  are  given  prior  to  discontinuing  the  run. 


• ILLEGAL  RESTRAINT  ON  SLAVE  NODE  --  This  means  a 

slave  node  has  restraints  in  the  Z or  0 directions  or  it  is  re- 
strained in  the  R direction,  but  the  radial  coordinate  is  not 
equal  to  the  radial  coordinate  of  the  first  or  last  master  node. 


This  also  occurs  if  node  is  designated  as  a slave  node  more  than 
one  time . 


• ILLEGAL  SLAVE  ELEMENT  --  A slave  element  does  not 
have  exactly  two  slave  nodes. 

• ILLEGAL  REQUEST  FOR  NODE  SHAPE  --  A request  is 
made  for  a special  node  shape  which  does  not  have  an  iden- 
tification number  of  1,  2,  3,  4 or  5. 

• ILLEGAL  REQUEST  FOR  ELEMENT  SHAPE  --  A request  is 
made  for  a special  element  shape  which  does  not  have  an 
identification  number  of  1,  2,  3,  4 or  5. 

• MINIMUM  TIME  INCREMENT  HAS  BEEN  VIOLATED  -- 
This  means  the  time  increment  given  by  Equation  (43)  is  less 
than  specified  by  DTMIN  on  the  Integration  Time  Increment 


. 


card.  The  results  are  written  on  Tape  2 and/or  Tape  3 if 
previously  specified  (they  are  not  printed).  If  it  is  desired 
to  continue  the  run  to  later  times,  DTMIN  must  be  decreased 
for  the  next  restart  run. 

• CP  TIME  LIMIT  EXCEEDED  — This  means  the  central  pro- 
cessor time  is  greater  than  specified  by  CPMAX  on  the 
Identification  card.  The  results  are  then  printed  and  written 
on  Tape  2 and/or  Tape  3 if  previously  specified. 

• ENERGY  CHECK  INDICATES  NUMERICAL  INSTABILITY  — 

This  means  the  kinetic  energy  of  the  system  is  at  least  10 
percent  greater  than  the  initial  energy  (kinetic  plus  internal 
explosive  energy)  of  the  system.  Since  this  is  not  physi- 
cally possible,  the  run  must  be  numerically  unstable.  The 
results  are  printed  and  written  on  Tape  2 and/or  Tape  3 if 
previously  specified. 

• REQUESTED  CYCLE  NOT  FOUND  ON  TAPE  --  This  means 
a restart  tape  (Tape  2)  does  not  have  data  for  cycle  number 
specified  on  the  Identification  card. 

4.  7 CENTRAL  PROCESSOR  TIME  ESTIMATES 

For  various  problems  run  on  a Honeywell  6080  computer,  the  EPIC-2  pro- 
gram has  used  approximately  0.  0045  central  processor  second  per  node 
per  cycle. 
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CENTRAL  MEMORY  STORAGE 
REQUIREMENTS  AND  ALTERATIONS 


The  EPIC-2  program  is  currently  limited  to  10  materials  (five  solid-liquid 
and  five  explosive),  1000  nodes,  2000  elements  and  five  sliding  surfaces, 
each  with  a maximum  of  100  master  nodes,  100  slave  nodes  and  100  slave 
elements.  This  requires  76K  Central  Memory  Storage  on  a Honeywell 
6080  computer.  These  limits  can  be  altered  by  altering  the  appropriate 
common  arrays  in  all  subroutines  which  contain  those  arrays. 
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APPENDIX  A 
SAMPLE  PROBLEM 


i 


Input  data  tor  a sample  problem  is  provided  in  Figure  A-l,  and  a geometry 
plot  is  shown  in  Figure  A-2.  This  problem  consists  of  a copper  rod  with  a 
hemispherical  nose  impacting  a steel  plate  at  HO,  000  inch/sec.  The  total 
length  of  the  rod  is  4.0  inches  and  the  radius  is  0.0  inch;  the  thickness  of  the 
plate  is  0.0  inch  and  the  radius  is  4.0  inches.  The  nodes  in  the  rod  are  gen- 
erated in  two  sections  to  allow  expanded  spacing  near  the  aft  end  of  the  rod; 
likewise,  the  nodes  in  the  plate  are  generated  in  two  sections.  The  nodes  in 

the  forward  end  of  the  rod  are  designated  as  slave  nodes,  and  the  nodes  on 

i 

the  top  surface  of  the  plate  are  designated  as  master  nodes.  This  problem 
contains  072  nodes  and  000  elements. 
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Figure  A-l.  Input  Data  for  the  Sample  Problem 
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